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FOREWORD 

This document is Volume 1 of a three volume report describing the 
Reacting and Multi- Phase (RAMP2) computer code developed by the Computa- 
tional Mechanics Section of Lockheed's Huntsville Research & Engineering 
Center, Huntsville, Alabama. Volume I describes the overall contractural 
effort and the theory and numerical solution for the computer code. Volume 
II provides a detailed description of all the elements used in the RAMP2 
code, ana Volume III is the program user's and applications manual. 

Documentation of the computer code was prepared in partial fulfillment 
of 'Contract NAS9-16256 with the NASA-Lyndon B. Johnson Space Center, 

Houston, Texas. The contracting officer's technical representative for this 
study was Mr. Barney B. Roberts, ET41. 

The author ackuowledges the efforts of Dr. Terry F. Greenwood of NASA- 
Marshall Space Flight Center and Mr. S.J. Robertson of Lockheed-Huntsville , 
both of whom contributed to the development of the RAMP2 code. 

Companion documents to this report Include a users and applications 
manual for RAMP2 computer code; a computer program maintenance manual for 
RAMP2; a report which describes "he modifications made to the NASA-Lewis 
TRAN72 computer code; the original documentation of the NASA-Lewis TRAN 72 
computer code; and the origins) documentation of the Bourdary Layer Integral 
Matrix (BLIMrJ) computer code. These documents are, respectively: 

• "High Altitude Chemically Reacting Gas-Particle Mixtures - Volume II 
- Program Manual for RAMP2," LMSC-HREC TR D86 7400-1 I. 
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• "High Altitude Chemically Reacting Gas-Particle Mixtures - Volume 
111 - Computer Code User's and Applications Manual, " LMSC-HREC TR 
D867400-III. 

• "User's Guide for TRAN72 Computer Code Modified for Use with RAMP 
and VOFMOC Plowfield Codes," LMSC-HREC TM D390409. 

• Svehla, R.A. , and B.J. McBride, "FORTRAN IV Computer Program for 
Calculation of Thermodynamics and Transport Properties of Complex 
Chemical Systems," NASA TN D-7056, January 1976. 

• Evans , R.M. , "Boundary Layer Integral Matrix Procedure BLIMP-J 
User's Manual," Aerotherm UM-75-64, July 1975. 
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Symbol 
A i 


a 



D/Dt 

dv 

A 

e 

E 

F 


f 


G 

H 

h 

I 

sp 

K 

k 


SYMBOLS AND NOTATION 
Description 

2 

parameter defined by Eq. (3.46) ,1/sec; area, ft 

speed of sound, ft/sec 

parameter defined by Eq. (3.86c) 

parameter defined by Eq. (3.83), ft /sec /°R 

drag coefficient, dimensionless 

specific heat at constant pressure, ft^/sec^/°R 

specific heat at constant volume, ft /sec /°R 

drag force, lbf 

substantial derivative 

3 

an element of volume, ft 
unit vector, dimensionless 

energy in the gas -particle system control volume, 

ft z /sec2 

interpolation factor, dimensionless; or total force 
acting on the system defined by Eq. ( 3 . 30) t lbf 

drag coefficient parameter (C~/Cjj ) , 

dimensionless otokes 

Nusselt number parameter defined by Eq. (3.77) 
dimensionless 

2 2 

total enthalpy, ft /sec 

7 7 

local enthalpy, ft /sec 
specific impulse, lbf-sec/lbm 

particle heat transfer film coefficient defined by 
Eq. (3. 76), lbm/sec/°R/ft 

thermal conductivity of gas, lbm/sec/°R 
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NG 


NP 


NS 


Nu 

n 

P 

P. 

1 

Pr 

Q 


q 

R 

Re 

rJ 

S 


T 

t 

u, v 

V 

W 



index denoting number of discrete species 
considered = NG-tNP 

Nusselt number, dimensionless 

unit normal vector 
2 

pressure, lbf/ft 

2 

partial pressure of species i, lbf/ft 
Prandtl number, dimensionless 

heat transferred to or from a gas -particle system 
control volume, Btu 

velocity, ft/sec 

gas constant, ft^/sec^/°R 

Reynolds number, dimensionless 

particle radius , ft 

surface area of a gas -particle system control 
volume; ft 2; entropy, ft^/sec^/oR 

temperature, °R 

time, sec 

velocity components, ft/sec 

3 

volume of a gas -particle system control volume, ft 

work performed on or by the gas -particle system 
control volume, ft^/sec2 
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Symbol 

w 

x, r 
Greek 
a 

P 

y 

6 

6 

T 

e 

V 

0 

M 

v 

X 

P 

a 

a 

L 

X 

8 


Description 

production rate, Ibm/sec 
position coordinates, ft 


accommodation coefficient, dimensionless; or Mach 
angle, radians 

characteristic slope, radians 
specific heat ratio, dimensionless 
any extensive quantity of the element dv 
Kronecker delta 

0, i for two-dimensional or axisymmetric flow, respectively 

2 

viscous stress tensor, lbf/ft 
emissivity, dimensionless 
nabla 

inclination of the flow vector with respect to the 
x-axis, radians 

chemical potential, cal/gm 

2 

viscosity, lbf-sec/ft 
equation modifier 

3 

density, slug/ft 

2 

surface stress tensor, lbf/ft 

2 3 

Stefan -Boltzman constant, ft /sec 
particle stream function, lbm-sec/ft 


indicates summation, 



j = l 

species mass fraction 
indicates partial derivative 
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Subscripts 

I 

i 

L 

m 

n 

w 

*. y 


Superscripts 


j 


Description 

indicates initial data surface 

indicates the quantity pertaining to species i 

denotes local surface conditions 

indicates the quantity pertaining to the gas -particle 
mixture 

indicates n*^ data surface 
denotes nozzle wall conditions 

denotes partial differentiation in the x and y (or r) 
directions 

denotes a vector quantity 

denotes an average value over a step length 
indicates the quantity pertaining to a particle specie 
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1. INTRODUCTION AND SUMMARY 

The overall objective of this study was to deliver a nozzle-plume flow- 
field code that has capabilities which do not presently exist in a single 
computer code. The RAMP code (Refs. 1 and 2), developed by Lockheed under 
government funding, was chosen as the basic code from which to work. The 
basic RAMP code is of modular construction and has the following 
capabilities: 

• Two-phase with two-phase transonic solution 

• Two-phase, reacting gas (chemical equilibrium reaction kinetics), 
supersonic inviscid nozzle/plume solution 

• Operational for inviscid solutions at both high and low altitudes. 

During the course of the study the following capabilities were added to the 
code to produce the RAMP2 program: 

• Direct interface with JANNAF SPF code (Ref. 3) 

• Shock capturing finite difference numerical operator 

• Two-phase, equilibrium/frczen, boundary layer analysis 

• Variable oxidizer- to-fuel ratio transonic solution 

• Improved two-phase transonic solution 

• Two-phase real gas semi-empirical nozzle boundary layer 
expansion 

• Continuum limit criteria 

• Sudden freeze free molecular calculation beyond the 
continuum limit 

• Interface between the Lockheed Plume Impingement Program 
(PLIMP) Ref. 4. 
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Most of the above capabilities already existed in other computer codes. 

These codes were incorporated into the RAMP code to enhance its usefulness. 
The emphasis of this effort was to improve the capability of the program to 
treat phenomena which significantly affect high altitude plumes. 

The contents of this volume along with the other two volumes (Refs. 4 
and 5) are designed to enable the user to easily, efficiently, and accu- 
rately apply the RAMP2 codes to calculate rocket exhaust nozzle and plume 
flow fields to solve practical engineering design problems. 

RAMP2 consists of three basic computational modules: the TPAN72 
program for generating equilibrium thermodynamic and transport data, the 
RAMP2F code for solving the flow field, and the BLIMPJ code used to 
calculate the nozzle boundary layer. Because of computer storage 
limitations the three programs must be executed separately; however, 
communication between the programs has been provided via temporary files so 
that, except for separate executions, they can be considered as one program. 

Detailed descriptions of the TRAN72 (Ref. 6) and BLIMPJ (Ref. 7) have 
not been included, although Volume II contains a brief description of the 
subroutines of each of the codes. Section 4 of Volume III contains informa- 
tion on how to use and input the TRAN72 code, and Section 5 contains a brief 
discussion of the BLIMPJ code. Complete descriptions of the two codes are 
available in Refs. 6 and 7. 

Solid propellant motors are frequently utilized to provide .aunch 
vehicle boost and stage separation. Most solid rocket motor propellants 
contain metal additives which increase the energy content of the system and 
suppress combustion pressure instabilities. The presence of metal additives 
results in condensed products in the nozzle-exhaust flow fields. This can 
create a number of adverse effects. 
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Condensed products (i.e. f particulate aluminum oxide, etc.) being inert 
can do no expansion weak. Consequently thr particles are accelerated 
through the flow field via drag exerted on the gas by the viscous shearing 
action between the gas and particulate phases. Since the particles do no 
expansion work, the gas phase cools more rapidly than the particles. The 
particles at any given location are thus at a higher temperature than the 
gas so that heat is thus transferred from the particles to the gas by 
conduction and radiation. The net result is that the gas phase expands 
useful work in accelerating the particles while acquiring heat from the 
particles. This is an irreversible, non-adiabatic process in the gas phase 
expansion. 

Exhaust plumes create hostile environments to surfaces immersed in the 
flow. These are in the form of structural loads, contamination, and 
heating. Heating results from both gas and condensates impinging on the 
surfaces. Particulates in the flow can also cause surface erosion and other 
structural damage. 

Exhaust plume applications and plume impingement problems typically 
occur during launch, staging, and rendezvous. The current Space Shuttle 
design offers several illustrations of areas where exhaust plume flowfield 
properties and flow structure must be known. This is schematically illus- 
trated in Fig. 1-1 which indicates potential plume related problems from 
solid rocket motors. During the Space Shuttle launch the solid rocket 
booster motor exhaust plumes impinge on the launch stand hardware. While 
the launch vehicle is in the vicinity of the launch pad, reradiation from 
the launch pad to the Orbiter i3 a potential problem. The solid rocket 
motor exhaust plume affects the Shuttle base drag and aerodynamics during 
the boost phase. The solid rocket separation motors mounted on the booster 
motor’s forward and aft ends that are used to effect stage separation can 
potentially subject the Orbiter and external hydrogen/oxygen tank to a 
rather severe environment. 
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Mission Profile Illustrating Rocket Plumes Associated with Various 
Shuttle Operations 
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Exhaust Impingement applications require a detailed definition of the 
gaseous plume structure as well as particle trajectories and dynamic prop- 
erties. Plume simulation studies require that plume shapes be known. The 
exit plane pressure must be known along with the gas thermodynamics. To 
adequately describe the flowfield properties, the mutual effect of gas on 
particle and particle on gas must be calculated. Also "real gas” effects 
can be significant and should be Included in the gasdynaic considerations. 

Solutions for the supersonic flow of a gas-particle mixture follow one 
of two approaches: (a) a fully coupled solution in which momentum and 

energy are exchanged between the gas and particle phases (Refs. 8, 9, and 
10), and (b) an uncoupled solution (Ref. 11) in which particle trajectories 
are traced through nozzle-exhaust plume flows. The uncoupled solution 
considers ”real gas" effects; however, it treats only gas effects on the 
particle and is more applicable for low aluminum content propellants. The 
coupled solutions are primarily performance oriented. These solutions 
readily permit treatment of highly aluminized propellants but are restruc- 
tured to constant thermodynamic properties. The performance code developed 
by Kllegel (Ref. 10) was subsequently utilized to trace liquid droplets in 
predicting contamination to surfaces (Ref. 12). 

Applications discussed previously require a knowledge of the flowfield 
structure. Previous studies (Refs. 13, 14, and 15) indicated the need to 
Include the treatment of "real gas” effects in nozzle-plume calculations. 
The decision was subsequently made to extend the coupled gas-particle 
solution (Refs. 8 and 9) to include the treatment of gas thermochemistry 
(chemical equilibrium and chemical kinetics). The choice of computation 
scheme and computer code becomes more obvious after examination of antici- 
pated applications. 
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The metnod-of-characteristics has proven to be a reliable numerical 
solution for nozzle-plume applications (Ref. 16). Existing gas-particle 
performance codes utilize the method of characteristics which indicated that 
an operational nozzle-exhaust plume code could be obtained with a minimum of 
numerical solution development. The gas-particle capability has been in- 
corporated into a streamline-normal code (Ref. 17) which utilizes the method 
of characteristics to solve for local flow properties. There are several 
reasons for choosing this approach; 

• Numerical difficulties encountered in highly expanded 
flow are circumvented. 

• The streamline-normal technique has a built-in 
mechanism for tracing particle streamlines. 

• Transition flow between the continuum and free 
molecular flow regime is more readily treated. 

The code, in its present form, will handle the flow problems which 
exhibit or have any of the following characteristics: 

• Supersonic Inviscid Flow 

• Highly Underexpanded Nozzles 

• Highly Overexpanded Nozzles 

• Shock Waves 

• Sliplines 

• Solid Walls 

• Pressure Boundary 

• Nonequilibrium Chemistry (Finite Rate) 

• Equilibrium Chemistry 

• Ideal Chemistry 

• Free Molecular Flow 

• Two-Phase Flow 

• Oxidizer-to-Fuel (0/F) Ratio Gradients 

• Nozzle Boundary Layer. 
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This report summarizes the improvements which were made under this 
contract to increase the capability of the RAMP program. In addition, this 
report presents a detailed development of the equations governing the super- 
sonic flow of a chemically reacting gas-particle mixture. Development of 
the governing relations follows to a large extent the work of Kliegel. 

Basic assumptions are stated and the set of partial differential equations 
describing the gas-particle system are developed. These relations are then 
cast in "characteristic” form and the corresponding difference equations 
written. Details of the numerical solution are described for the various 
data point types. The presentation is then concluded with a description of 
the numerical integration of the conservation equations. 

In addition to a description of the above analysis techniques, appen- 
dixes are Included which discuss: (1) particle drag and heat transfer 

coefficients; (2) non-isoenergetic gas-phase flow treatment; (3) chemical 
equilibrium calculations in gas-particle flows; (4) non-continuum flow 
expansions; and (5) integration of the finite-rate chemical kinetic 
equations. 
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2. SUMMARY OF STUDY PERFORMED FOR CONTRACT NAS9-16256 

A precise knowledge of local flow properties in nozzles and exhaust 
pluses is necessary for performance, radiation, attenuation, heat transfer 
and isplngeoent analyses. The reacting and multiphase (RAMP) computer 
program is designed to give detailed flowfield information in the supersonic 
region of a reacting multi -phase, two-dimensional or axisymotetric flow 
field. The boundaries of the flow field may be solid such as in a nozzle or 
"free" such as in a plume. The analysis may be utilized therefore to pre- 
dict performances as well as plume characteristics of a given engine system. 
A printed record of the program results is given for user inspection while a 
binary tape is provided for subsequent manipulation by other analyses. 

The flow of a gas-particle mixture is described by the equations for 
conservation of mass, conservation of momentum, and conservation of energy. 
In the gaseous phase the state variables P, p, R and T are related by "he 
equation of state while for the particulate phase the equations are for the 
particle drag, particle heat balance and the particle equatii of state. 
Development of these equations is based on the following assumptions: 

1. The particles are spherical in shape. 

2. The particle internal temperature is uniform. 

3. The gas and particles exchange thermal energy by 
convection and radiation (optional). 

4. The gas obeys the perfect gas law and is either frozen 
and/or in chemical equilibrium, or is in chemical 
non-equilibrium . 

5. The pressure of the gas and the drag of the particles 
contribute to the force acting on the control volume. 
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(* . The gas Is inviscid except for the drag it exerts on the 
particles. 

7. There are no particle interactions. 

3. The volume occupied by the particles is negligible. 

9. There is no mass exchange between the phases. 

10. A discrete number of particles, each of different size 
or chemical species is chosen to represent the actual 
continuous particle distribution. 

11. The particles are inert. 

The supersonic two-phase solution accepts the starting line provided by 
the internally calculated transonic solution as well as other pertinent data 
supplied through the read function. The equations of motion under the 
assumptions just listed are hyperbolic and permit the use of a forward 
marching scheme; a streamline/normal grid structure is employed where the 
step lengths in the axial and radial directions are under program control. 
Both BCD (printer) and unformatted binary output tapes are produced. A 
Frandtl-Meyer expansion of the gas phase and a free boundary calculation are 
employed to treat the plume flow solution. The run is terminated when pre- 
tpecifled problem limits are reached. 

The two-phase flow analysis will treat an extremely wide range of 
operating conditions. With few exceptions the limitations are imposed by 
the theory rather than numerical considerations. In this discussion 
dimension statement sizes which are arbitrarily set are not considered a 
limitation. Th'_ true limitations are: 

• S* ersonic regions influenced by embedded subsonic regions. 

Vacuum or limiting expansion limitation - a small region of the 
expansion fan for a vacuum expansion cannot be treated where the 
Mach number is so large that treatment by continuous flow assump- 
tions becomes meaningless (this limitation is both numerical and 
theoretical) . 

• Fo* two-phase flow the lower boundary can only be horizontal (l.e., 
..3zzle centerline). 


2-2 


LOCKHEED- HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC TR D867400-I 


During the course of the study the following capabilities have been 
added to the code. 

• Direct interface with JANNAF SPF code 

• Shock capturing finite difference numerical operator 

• Two-phase, equilibrium/frozen, boundary layer analysis 

• Variable oxidizer-to-fuel ratio transonic solution 

• Improved two-phase transonic solution 

• TWo-phase real gas semi-empirical nozzle boundary layer 
expansion 

• Continuum limit criteria 

• Sudden freeze free molecular calculation beyond the 
continuum limit. 

• Interface with Plume Impingement models, and 

• Documentation. 

2.1 MODIFICATIONS TO THE RAMP CODE 

This section deals with the modifications which were incorporated into 
the RAMP code under this contract. 

2.1.1 JANNAF Standard Plume Flowfield (SPF) Code Interface 

To perform a plume calculation, the SPF code (Ref. 3) requires nozzle 
exit properties as Initial conditions. The RAMP code was modified to punch 
or put on tape exit plane data in the format that the SPF code uses. The 
code will punch data for both single and two-phase cases. 

At present, the SPF code can accept a single chemical species distribu- 
tion which must be applied at all points across the exit. The RAMP code can 
have a distribution of species, so it was necessary to mass flow average the 
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species distributions across the exit plane to produce a single set of 
species. This restriction can be easily removed whenever SPF haa the capa- 
bility to handle spatial variations in chemical species. 

The code is set up to punch the species information for the various 
chemical systems that are available In the SPF code. The user must specify 
which system to punch. If the thermodynamic data are input to the RAMP code 
using cards then the species mole fractions that are to be punched for the 
SPF code must be input for the appropriate system in the same order as 
listed in Table 2-1. For finite rate cases these species are determined 
inside the code and need not be input. The user may also input another 
system. Additionally, ideal gas cases can also be handled. 


Table 2-1 

SPF CHEMICAL SYSTEMS 


System 1 — Hydrogen/Oxygen 

h,h 2 ,h 2 o,n 2 ,o,oh,o 2 

System 2 — Carbon/Hydrogen/Oxygen 

co, co 2 , h, h 2 , h 2 o, n 2 , o, oh, o 2 

System 3 — Carbon/Hydrogen/Oxygen/Chlorine 

A! 2 0 3 (S), CO, C0 2 , Cf, Cl 2 , H, H 2 , H 2 O t HCf, N 2< O, OH, C 2 

System 4 — Carbon/Hydrogen/Oxygen/Chlorine/Fluorine 

CO, C0 2 , Cf, Cf 2 , F, F 2 , H, H 2 , H z O, HCf, HF, N 2 » O, OH, 0 £ 

System 5 — Hydrogen/Oxygen/Boron 

BO, B0 2 , B 2 0 3 , H, H 2 , H 2 0, HB0 2 , N 2 , 0, OH, 0 2 

System 6 — Hydrogen/Oxygen/Boron/Chlorine/Fluorine 

BF, BF 2 , BF 3 , BO, BOCf, BOF, B0 2 , B^, CO, CO z , Cf, Cf 2 
F, F 2 , H, H 2 , H 2 0, HB0 2 , HCf, HF, N 2 , O, OH, O z 
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2.1.2 Shock Capturing Finite Difference Operator 

The original RAMP code has the logic for computing any number of right 
or left-running shocks using shock fitting techniques. This logic was 
partially checked out under previous efforts. It is desirable to have the 
capability to treat shocks that can occur in some nozzles. Shock capturing 
schemes require no special equations or program logic and are reliable. For 
most nozzle flows, existing shock capturing techniques are sufficiently 
powerful to treat the shocks. Additionally, in order to directly interface 
with the SPF code, it would be desirable to use a shock capturing numerical 
operator. 

For the above reasons a shock capturing algorithm was added to the 
code. The methodology, equations and grid system which was incorporated 
into the code is identical to the SPF code (Ref. 3). To obtain more inform- 
ation on the scheme the reader is referred to Ref. 3. 

For the shock capturing option, the computational grid is composed of 
vertical lines and a radial coordinate system normalized by the local wall 
coordinate. This approach as taken to alleviate the problem of inserting a 
point per line for a Cartesian system. The computational procedure is as 
follows: the forward marching algorithm is used to compute the flowfield 
properties at all interior points. Flowfield properties at the wall are 
obtained via the method-of-characterlstics solution. Once the wall has been 
computed, the coordinate system is transformed and the axis and interior 
points are solved. 

To check out the validity of the shock capturing option, comparisons 
were made with known method-of-characteristlcs solutions. Figure 2-1 
presents pressure distribution along a nozzle centerline and wall for 
various calculatlonal techniques. As can be seen from the figure, the shock 
capturing solution is identical to the characteristics and SPF solutions. 


2-5 


LOCKHEED-HUNTSVILLE RESEARCH & ENGINEERING CENTER 



Mach Number 


LMSC-HREC TR D867400-I 



0 5 10 15 20 25 


Axial Distance from Nozzle Throat (Throat Radii, 


Fig. 2-1 Nozzle Wall and Centerline Mach Number Distributions 
Using Various Calculation Techniques 
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2.1.3 IVo-Phase, Equilibrium/ Frozen, Boundary Layer Analysis 

Nozzle boundary layers are known to influence certain regions of nozzle 
flow and high altitude exhaust plumes (Refs. 18, 19 and 20). For nozzle 
designers, the nozzle boundary layer is important in determining the thermal 
loads to the nozzle, performance losses due to heat transfer and effects on 
the nozzle pressure distribution due to the displacement thickness effect on 
the inviscid flow structure. For spacecraft designers the nozzle boundary 
layer is important because of its effect on the exhaust plume. At high 
altitudes the nozzle boundary layer causes the plume to expand to large 
angles (approximately 180 deg). In these backflow regions, spacecraft and 
sensitive surfaces are subjected to unwanted contamination, forces, moments 
and heating rates. For some applications the radiative properties of the 
expanded boundary layer flow is importac. . Thus, it is easy to see that for 
many applications the nozzle wall boundary layer is an important factor. 

To date there have been numerous techniques used to account for 
boundary layer effects of the nozzle and exhaust plume flow fields. 
Lockheed-Huntsville has an option in the MOC and PIP code (Ref. 21, a 
version of RAMP) which superimposes a power law boundary layer profile based 
on a flat plate boundary layer thickness on the exit plane start line. The 
profile can contain a total temperature variation which has been found to be 
Important (Ref. 21) in the back flow region of the plume. This is a very 
simplified approach and the method of determining the boundary layer thick- 
ness is very crude. Lockheed-Huntsville personnel and other Investigators 
have used various boundary layer codes and the inviscid nozzle results to 
generate a better set of boundary layer properties which were then super- 
imposed on the nozzle exit plane Inviscid results. Next, some investigators 
have gone a step further and iteratively solved the inviscid nozzle and 
boundary layer to determine the effect of the boundary layer on me nozzle 
flow field. Finally, it is possible to do a fully viscous nozzle solution 
with a method such as Lockheed-Huntsville's GIM code (Ref. 22). However, a 
fully viscous nozzle solution is beyond the scope of this effort. It should 


2-7 


LOCKHEED-HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC TR D867400-I 


be noted that all techniques except the first require a tremendous amount of 
effort and hand calculations to finally obtain a plume solution. The pur- 
pose of this task was to automate the boundary layer nozzle solution to 
arrive at an exit plane start line which Includes all the effects mentioned 
above. 

The BLIMPJ boundary layer code (Ref. 7) was chosen for the solution of 
the nozzle wall boundary layer. BLIMPJ is a JANNAF standard boundary layer 
code for determining boundary layer effects on the performance of a rocket 
engine. The BLIMPJ code is set up to handle nozzle flows, equilibrium or 
frozen chemistry, uses the JANNAF standard thermochemistry curve fits and 
has numerous ways of treating the nozzle wall that allow it to treat liquid 
engines as well as solid motors with and without ablative walls. At 
present, the Interface between RAMP and BLIMPJ allows only the following 
wall options: adiabatic wall, steady state energy balance or a specified 

temperature distribution. The adiabatic wall option is the default option 
in the code. The remaining wall options in the BLIMPJ code could be used 
with a few modifications to the RAMP code. In addition to wall boundary 
condition options, the user may also select laminar, turbulent, frozen, or 
equilibrium flows or let the code default to turbulent, equilibrium. If the 
turbulent boundary layer option is used, the flow will transition to a 
turbulent profile when the Reynolds number based on momentum thickness 
exceeds 250. 

All input data for the boundary layer solution are generated inside the 
RAMP code except the nozzle wall temperature distribution, if that option is 
specified. Chemical species must be input if the thermochemistry data for 
the inviscid solution was input via cards. 

The BLIMPJ code was too large for conversion to a module for the RAMP 
code, so if it is desired to generate a nozzle boundary layer and obtain a 
viscous exit plane start line it is first necessary to run the RAMP code to 
calculate the inviscid nozzle solution. The RAMP code generates a file with 
the BLIMPJ input data on it. The BLIMPJ code is then executed (no data 
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cards). The BLIMP J results are printed out and stored on a file for sub- 
sequent use by the RAM? code. The RAMP code is then reexecuted with just a 
single input card and an exit plane vertical or normal start line is gene- 
rated. The exit plane start Tire has the results of the inviscid nozzle and 
boundary layer solutions merged. This capability is a considerable Improve- 
ment over previous methods which required many hand calculations, cross 
plotting, etc. 

For most rocket motors the boundary layer is fairly thin (approximately 
5 to 10% of exit diameter) and the resultant effect on invisci'’ flow proper- 
ties is minimal. For these cases one pass through the inviscid nozzle 
solution and boundary layer calculation is adequate. The boundary layer 
results will then be superimposed on the inviscid nozzle solution and an 
exit plane start line with boundary layer effects will be generated which 
can be used to perform a plume expansion. 

Some low pressure or low thrust motors have boundary layers which 
contain a significant portion of the total mass flow of the system. In 
these cases the entire solution should be iterated for by making two passes 
through the inviscid nozzle and boundary layer calculations. After the 
initial nozzle and boundary layer solution has been completed, the nozzle 
solutions will be calculated with the actual wall contour adjusted by the 
local displacement thickness which was determined from the boundary layer 
calculation. The boundary layer solution will then be rerun using the new 
edge conditions from the second nozzle calculation. The results will then 
be superimposed on the original nozzle contour and second flowfield solution 
and an exit startline will be output or saved. This capability is not 
presently operational. However, for these cases the boundary layer solution 
is solved a second time with new edge conditions taken from the inviscid 
results at the edge of the first BLIMPJ boundary layer. Thus, the inviscid 
and viscous results match very well at the boundary layer edge. This option 
is controlled Internal to the program and requires no user interface. 
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2.1.4 Particle Tracing Through Boundary Layer 

Any par«.i.c’ late matter (A^O^-solids or unburned propellant drop- 
lets In liquid motors) which might enter the nozzle wall boundary layer 
could be influenced by the boundary layer to be expanded into the backflow 
region. Thus, there may be certain motor/nozzle configurations in which it 
is necessary to track particles through the boundary layer in order to 
obtain the most accurate representation of the nozzle flow properties at the 
exit plane. 

Lockheed-Huntsvllle added a particle streamline tracing module ' the 
nozzle code. The code which was used as the basic building block or this 
module was used in predicting the 1US, SSUS particle distributions published 
in Ref. 23. This code uses initial particle properties (velocity, flow 
angle, temperature and density) and traces the particles through a known 
flow field. This option is user-selectable. The basic method used is shown 
in Fig. 2-2. 

After the invlscld nozzle solution and boundary layer solutions hav* 
been completed, the RAMP code is re-executed and logic is banded off to the 
particle trajectory tracing module. The first step is to search the 
boundary layer edge location in the inviscid nozzle solution to determine if 
any pan icles are found to penetrate the boundary layer. Then for each 
particular Blze group which has been found to penetrate the boundary layer a 
distribution of individual particle properties will be determined (inter- 
polated off nozzle solution tape) at several stations along the boundary 
layer edge. This will be done for each size group. Next a map of the 
boundary layer properties will be generated using the boundary layer data 
tape. Finally, each particle streamline will be traced through this mapped 
flow field and particle properties will be determined at the previously 
determined exit plane. The effect of the particles on the boundary layer 
flow properties due to energy and momentum interchange are not considered. 
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Particle 

Streamlines 



1 


Fig. 2-2 Schematic ,<f Particle Trajectory Tracing Method to be Used 
In the Invlscld Nozzle Code 
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Finally, the boundary layer results, the inviscld nozzle solution and 
the particle properties in the boundary layer are merged- and an exit plane 
start line for a RAMP restart or SPF run will be generated. 

Several cases have been executed using the particle trajectory option. 
Figure 2-3 presents some results of the Large Boeing IUS motor calculations 
which illustrate the effect the boundary layer can have on the particle 
properties. The axial distribution of a particle limiting streamline 
temperature and velocity are shown for both the inviscid and viscous 
calculations. 

2.1.5 Boundary Layer Expansion at Lip 

Another required capability of the nozzle code is a fully automated 
technique to expand the boundary layer at the lip. After an exit plane 
start line has been generated it is very straightforward to use the Prandtl- 
Meyer relationship to expand the flow. The key point here is to get a fully 
supersonic start line. Numerous investigators have made various assumptions 
to obtain start lines for plume solutions that account for the nozzle wall 
boundary layer. The various techniques will be discussed briefly. 

Perhaps the most simple method of obtaining a supersonic start line is 
to discard the subsonic portion of the boundary layer (Refs. 18 and 24) and 
adjust the spatial distribution of boundary layer points up to the lip. 

This method throws mass away from the system and will affect the results in 
the backflow region. For nozzles which have thin turbulent boundary layers, 
these effects are probably minimal. For very viscous nozzles (low pressure) 
or laminar boundary layers the effects of this assumption could severely 
compromise the plume results in the backflow region (greater than 30 deg). 

It should be noted that this technique could consider total temperature 
variations through the boundary layer. 
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Fig. 2-3 Large Boeing Motor 0.5 micron Radius Particle Limiting Streamline 
Temperature and Velocity distributions for Calculations With and 
Without the Nozzle Boundary Layer 
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Seubold (Ref. 25) used the technique of replacing the subsonic portion 
of the boundary layer vith mass average properties at a slightly supersonic 
value. This resulted in a layer of constant properties near the vail. The 
resultant plume flow had the correct mass flow but did not consider the 
total temperature variations. 

The first two techniques assumed that the static pressure through the 
boundary layer is constant. Near the lip, however, the subsonic portion of 
the boundary layer is influenced by the expansion process (Refs. 26, 27, and 
28) as shown in Fig. 2-4. For highly underexpanded flows the sonic line has 
been found to attach to the lip (Ref. 26) so that the flow at the wall must 
rapidly accelerate when it gets near the lip and the static pressure rapidly 
decreases. For overexpanded flows the subsonic flow merely stays subsonic 
as it negotiates the lip so that downstream flow conditions can feed back up 
into the boundary layer. 



Fig. 2-4 Mach Number Contours for an Underexpanded Axially 
Symmetric Nozzle (Case 6) 
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Cooper (Ref. 29) went a little further than the first two approaches as 
far as treatment of the subsonic flow expanding to the sonic li.> conditions. 
He assumed that the flow (M”0 at wall) just upstream of the lip expanded to 
M“1 at the lip. Since the static pressure upstreaa of the lip at the wall 
was the saae as the wall total pressure and the frees tream static pressure, 
the only way for M**l to be reached at the lip was for the wall static pres- 
sure to drop. Therefore, there will be a static pressure distribution 
across the flow at the lip. It was then assumed that the static pressure at 
a specified point in the boundary layer is not if fected by the expansion and 
the pressure distribution is faired into the known lip conditions. This 
aethod is not totally rigorous but does include boundary layer effects and 
has been shown to result in good comparisons with some experimental data in 
the backflow region. 

Finally, there are exact solutions to the corner flow problem, bird 
(Ref. 26) uses a direct simulation Monte Carlo aethod. He set his model up 
so that it would compute the entire nozzle or the region near the throat. 

He predicts the attachment of the sonic line for underexpanded nozzles. 

Baum (Ref. 28) uses a finite difference method along with boundary layer 
equations to describe the subsonic portion of the boundary layer. His 
application was for base flow about blunt bodies and compares well with 
data. Finally, the previously mentioned GIM code can be used to exactly 
solve the corner problem. All three of these methods ar?. very complex and 
are outside the scope of this effort. 

In the present version of the code the method of Seubold was incorpo- 
rated. Total temperature variations are included and the mass flow in the 
boundary layer is conserved. Plume calculations using this method were made 
for a 5 lbf bi propellant motor. The mass flux predictions as shown in Fig. 
2-5 compare very well with experimental data. 
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Fig. 2-5 Mass Flux Normalized by Centerline Mass Flux 
vs Angle from Plume Centerline for 5 lbf 
Bipropellant Motor 


2.1.6 Variable Oxldlzer-to-Fuel Ratio Transonic Solution 

Solution of the subsonic-transonic region of a liquid rocket engine can 
vary In complexity from a simple one-dimensional variable 0/F stream tube 
analysis (Ref. 30) to the most detailed model such as that of the Distri- 
buted Energy Release (DER, Ref. 31) model. The streamtube analysis performs 
a multizone one-dimensional calculation to the sonic point given a known 0/F 
distribution just downstream of the injection face. The DER program is a 
complex model which is initiated upstream of the injector face and continues 
the solution up through the sonic line. The DER code was used to initiate 
nozzle solutions in Ref. 32 but is not particularly easy to use and input 
and requires a good bit of user experience to successfully execute. For 
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these reasons it is not utilized in the nozzle code. On the other hand, 
a one-dimensional streamtube analysis does not account for two-dimensional 
effects. A time-dependent scheme (Ref. 33) is a compromise between these 
two schemes. The approach includes the radial momentum equation which 
results in a set of mixed partial differential equations. The solution 
procedure is an unsteady time-dependent finite difference technique with 
equilibrium chemistry. This technique has both the equilibrium and variable 
0/F chemistry option and has been utilized (Refs. 34 and 35) previously with 
excellent results. This code has been incorporated into the RAMP code as a 
module and has been executed for all combinations of ideal and equilibrium 
chemistry, constant and variable 0/F distributions. Additional inputs 
required by this module are: (1) combustion chamber geometry down to the 

throat; (2) average 0/F ratio; (3) location of the initial data surface to 
start the transonic solution; (4) location of the nozzle throat, and; (5) 
the 0/F distribution at the initial data surface. 

Figure 2-6 presents some results of the variable 0/F calculations for 
the Space Shuttle Reaction Control System engine. Exit plane distributions 



Radial Distance from Plume Centerline (in.) 


Fig. 2-6 Comparison of RAMP Calculation for Space Shuttle RCS 
Motor Plume with Pitot Pressure Survey Taken at 45 in. 
from Nozzle Exit Plane 
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of pressure and temperature are presented for both constant and variable 0/F 
assumptions* The constant 0/F calculation was performed assuming Mach 1*05 
at the throat. The variable 0/F calculation was initiated upstream of the 
nozzle throat at the point where the nozzle begins to contract. The 0/F 
distribution at this location was inferred from the injector pattern and 
incoming mass flow distribution, 

2.1.7 Improved Two-Phase Transonic Solution 

The original transonic module which was Incorporated into the RAMP code 
could handle throat radii of curvature to throat radius ratios above 1.5. 
Many solid motors have radii of curvature ratios smaller than 1. To 
alleviate this limitation the improved approximate transonic module was 
taken from the new Standard Performance Prediction program (SPP, Ref. 36) 
and put into the RAMP code. 

Many transonic solutions using this module have been run for varying 
upstream and downstream radii of curvature ratios. 

2.1.8 Continuum Limit Criteria 


There have been numerous studies over the past several years concerning 
methods of determining where continuum flow breaks down and free molecular 
flow begins. Examples of these criteria are the Knudsen number criteria, 
criteria based on Mach and Reynolds number, and the "breakdown parameter" as 
proposed by G.A. Bird (Ref. 37). 


The present version of the program uses an average Knudsen number to 
check for transition from continuum to "free molecular" flow regimes where 
the average Knudsen number between the old streamline base point is calcu- 
lated via the following equation: 


Kn » .788539 CM 2 /Rg) 


In T, 


In T 2 


/dS 
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where the properties are averaged between the old (1) and new (2) streamline 
points. The flow regime is determined by checking the calculated Knudsen 
number against the input Knudsen number for criteria translational freezing. 
Once the flow regime has been determined the appropriate specific heat ratio 
(gamma) and total conditions are calculated. A Knudsen number ot 10 is 
standardly used to check for transition. 

The criteria of Bird (Ref. 37) have also been included in the program. 
The Bird criteria specifies where the continuous flow equations start to 
break down. This is not necessarily the location in the flow where the flow 
goes free molecular. 

Bird's breakdown criteria 


P ,-i !£ 

V ds 

is calculated at each point on the normal. If the breakdown parameter of 
0.05 is exceeded, the code points out the location on the normal where P = 
0.05. The collision frequency v is calculated using the hard sphere 
relationship: 


5^ PRT 
V ' 4 U 

The code does keep track of the continuum "breakdown criteria" and 
locates the surface in the plume where the "breakdown criterion" of 0.05 is 
located. In the future the appropriate data at this surface will be stored 
so that this information may be passed along to a full Monte Carlo solution. 
An example of the application of the Knudsen number and Bird breakdown 
criteria is shown in Fig. 2-7. 
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Fig. 2-7 Large Boeing Motor Exhaust Plume Mach Number Contours 
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The Interim Upper Stage (IUS) is used to boost payloads from near earth 
orbit to higher orbits. Ibis stage uses solid rocket motors as propulsion. 
From a payload and stage design standpoint it is necessary to define the 
exhaust plume. The RAMP code is especially suited to calculating the IUS 
exhaust plume including all the driving phenomena: equilibrium/ frozen 

thermochemistry with total enthalpy variations, two-phase flow, nozzle wall 
boundary layer, and free molecular flow. Figure 2-7 presents contour plots 
for Mach number in the exhaust plume. Shown on this figure is the Knudsen 
number of 10 and the continuum breakdown surface corresponding to the Bird 
breakdown criteria of 0.05. Beyond the Knudsen number 10 contour the flow 
was treated as free molecular. 

2.1.9 Free Molecular Flow Option 

A "sudden freeze” free molecular flow option has been Incorporated into 
the code. Once the specified Knudsen number criteria have been exceeded the 
flow is frozen. (Specific heat ratio, molecular weight, velocity and 
temperature along <. gas streamline is held constant.) At this point the 
streamline is assumed to expand at a constant flow angle and the density 
varies according to the streamtube area. This option allows the calcula- 
tions to be performed in the backflow region (>90 deg) of the plume. 

Numerous cases have been run for several hundred nozzle exit diameters. 

2.1.10 Interface with Plume Impingement Model 

Specification of the amount and type of plume exhaust products at a 
surface Immersed in a rocket plume is Important for determination of space- 
craft surface degradation. After an adequate representation of the char- 
acteristics of the exhaust plume has been performed it is necessary to 
rela.' the spatial characteristics of the plume to a surface that is 
immersed in the plume. 
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Under previous studies, Lockheed-Huntsvllle developed the Lockheed 
Plume Impingement code (Ref. 38). This program will provide forces, moments 
and heating rates to bodies Immersed in the exhaust plume generated by the 
Lockheed Method-of-Characterlstlcs Program (Ref. 16). 

The Plume Impingement Program (PLIMP) has been modified so that it can 
use the results of the RAMP2 code to predict forces, moments and heating 
rates due to plume impingement. Additionally, the PLIMP code has been 
modified to determine the types and amounts of plume species at any given 
location of a body Immersed in the exhaust plume. This version of the PLIMP 
program is available with RAMP2. 
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3. FUNDAMENTAL EQUATIONS FOR STEADY FLOW 
OF REACTING GAS-PARTICLE MIXTURES 


3.1 BASIC ASSUMPTIONS FOR THE GOVERNING EQUATIONS 

The flow of a pas -particle mixture is described by the equations for 
conservation of mass, conservation of momentum and conservation of energy. 
In the gaseous phase the state variables, P,p, R and T are related by the 
equation of state while for the particulate phase the equations are for the 
particle drag, particle heat balance and the particle equation of state. De- 
velopment of these equations is based on the following assumptions: 

1. The particles are spherical in shape. 

2. The particle internal temperature is uniform. 

3. The gas and particles exchange thermal energy by con- 
vection and radiation (optional). 

4. The gas obeys the perfect gas law and is either frozen 
and/or in chemical equilibrium, or is .n chemical non- 
equilibrium. 

5. The forces acting on the control volume are the pressure 
of the gas and the drag of the particles. 

6. The gas is inviscid except for the drag it exerts on the 
particles. 

7. There are no particle interactions. 

8. The volume occupied by the particles is negligible. 

9. There is no mass exchange between the phases. 

10. A discrete number o f particles, each of different size or 
chemical species, is chosen to represent the actual con- 
tinuous particle distribution. 

11. The particles are inert. 

These are basically the same as originally stated by Kliegel except for the 
provision to: (a) calculate the radiation exchange between the gas and particle 
phase (assumption 3); (b) the gas phase can be either frozen, in chemical 
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equilibrium, or chemical non -equilibrium (assumption 4). The total mass 
and energy of the gas -particle system are not constant since provisions are 
made for particle streamlines to penetrate the exhaust plume boundary. 

3.2 GOVERNING EQUATIONS FOR THE GAS-PARTICLE MIXTURES 
3.2.1 Continuity Equation 

The various forms of the continuity equation for chemically reacting 
gas -particle mixtures are derived in this section. The forms of the equation 
derived, in order are: 

1. Species Continuity Equation (3.6) 

2. Global Continuity Equation /(Steady State) (3. 13)/(3.27) 

3. Gas Continuity Equation (3.19) 

4. Particle Continuity Equation (3.20) 

5. Species Continuity Equations as a Function of (3.25)/(3.29) 

Mass Fraction/(Steady State) 

6. Particle Continuity Equation for Each Particle (3.26)/(3.28) 

species j/(Steady State) 

Consider a chemically reacting gas -particle mixture flowing through 
some arbitrary stationary control volume, v, (Fig. 3-l)bounded by the con- 
trol surface, S. 



z 

Fig. 3-1 - Control Volume for the System of Equations 
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For such a flow system, Reynold's Transport theorem states that: 

£ / Ad « = / l? dV + /*<•“«; (3.1) 


V 


where A = any arbitrary parameter. 


Applying the Gauss, or Divergence theorem 

J B • ti dS = J V. Bd V ; 
S v 

where B = any arbitrary vector quantity 
to Eq. (3.1) yields: 


§t f A dV = / dV + / V . A? dV . 


(3.2) 


(3.3) 


Furthermore, the substantial or Euler's derivative may be written in the 
form: 

T5T = • VA • (3.0 

In a flow system of gas -particle mixtures in which r’.em cal reactions 
take place, the principle of conservation of mass of each chemical species 
may he written as: 


& /p 4 <iV - y* *. dV = 0 x = l, NS; 


(3.5) 


where w^ = production rate of species i due to either internal or external 
sources such as chemical reactions and mass additions. 


Applying Eq.(3.3)to the first term of Eq. (3.5) 

9P ; 


h / pi dV = /-9r dv+ / v -Pi«i dV ( = i.ns, 
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substituting the result back into Eq.(3.5) 



V • p, a. dV 
1 n 



V 


w. dV 
x 


0 i = 1, NS , 


and rearranging terms yields 




V • p^ q dV + 



i = 1, NS . 


The above equation simply states that, for species i, the mass rate of 
increase inside the control volume is equal to the net rate of mass flow into 
the control volume plus the rate of mass produced due to chemical reactions 
and mass addition. In the following, we shall, however, exclude mass addition 
from external sources. 


Combining the integrands under one integral 

+ V • P ; - w. 1 dV =0 i = 1 , NS , 


/ (lT + v-P i q i -w^ 


ana noting that the volume V under consideration is arbitrary; the only way 
for the above equation to be valid for all V is for the integrand to vanish. We 
therefore 1. ave 


9 P* 

—■ + v • P i - W. = 0 i = 1 , NS . (3.6) 

Equation (3.6) is known as the species continuity equation. It is valid 
for each chemical species at each internal quantum state. We shall, however, 
assume thal the various internal modes of motion are fully excited and are 
in equilibrium with each other. It is well known that this is approximately 
the case for the translational and rotational degrees of freedom where the 
equilibrium value is attained in a few collisions. In general, the vibrational 
degree of freedom aoproaches the equilibrium state somewhat more slowly, 
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except at very high temperatures. As chemical reactions usually occur at 
high temperatures, this approximation is often justified. The global continuity 
equation can now be derived by summing Eq. (3.6) for all species present. 


NS 


/QPi 

. \ 

2^(-ar + v.Pi q 

- w. = 0 

v 

i=l x 

/ 

NS NS 

NS 


-E** * 

i=l i=l 

i = l 


To arrive at the final form of the global continuity equation, it is necessary 
to take a closer look at the flow system. 


In the flow system analyzed, the time variation of the thermodynamic 
functions is slow compared to the longest relaxation time of the system. 
Therefore, the assumption that thermodynamic local equilibrium exists can 
be made. In such a system, the thermodynamic quantities for the nonequi- 
librium system are the same functions of the local state variables as the 
corresponding equilibrium thermodynamic quantities. Therefore, the partial 
specific quantity in an arbitrary volume element dV of a nonequilibrium 
system may be defined by the equilibrium relation 




(£) 


T, P, m. 


The quantities being held constant during the differentiation are the locally 
defined temperature T, pressure P, and the masses m. of the other NS-1 
species . 


The specific quantity <f> m (<$> per unit mass) is given in terms of the 
partial specific quantities by 


P 4> 

m 


NS 

X = ! , NS; 
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where the relationship between the mass density p of the mixture and the 

m 

species density p. is given by 

1 NS 

P m = ^ p. = H- id density (3. 


The specific velocity, q^, specific enthalpy, h^, and specific internal 


energy, are given by: 


NS 

^m ^i ^i 

i = l 

NS 

p h = V* p. h. 
m / „ *i i 

i = l 


p e r V 1 p. e. 

■m m / j l i 


Substituting Eqs.(3.8)and (3.s) into (3.7), 


+ V ‘ p m^m- = ° ’ 


and noting that for a closed system, the total mass for chemical reactions 


is conserved 


NS 

E»i * 

i = l 


Equation (3. 7) reduces to 


+ V • p q = 0 . 
'm mi 


Equation (3.13) is known as the global continuity equation. It can be written 
alternati vely as 

Dp m 

— ^m 7 - V = 0 

by means of the substantial derivative. 
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Equations (3,8), (3.9) and (3.12) can be rewritten in terms of gaseous 
and particle species as follows: 

Equation (3.8) 


NS 

f>m = E Pi 

i = 1 


Equation (3.9): 

NS 

^i ^i 
i = l 

Equation (3 . 12)*. 

NS 

E*i * 0 * 

i=l 


NP 

P + 23 ; (3.U) 

j=l 


NP 

- pq + 23 pj ^ j (3,15) 

j=i 


NP' 

v + ^ ^ w^=0 (3.16) 

j=l 


Assuming that reactions of the form A(gas) + B(gas)^z^ C(particle) + D(particle) 

9 

do not contribute significantly to the system (interchange of phases is negligible 
due to chemical reactions), Eq.(3.16) implies that 

N p 

w = 0 and - 0 . (3.17) 

j=l 


Substituting Eqs. (3. 14) through (3. 17) into Eq. (3.7) 


d_ 

9t 



+ V 



0 , 


and rearranging terms yields 


M 

dt 


+ v 



o 


(3.18) 
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Since the gaseous species and particle species do not interchange phases, 
both contributions must vanish. Therefore 


and 



-o 


(3.19) 

(3.20) 


Equation (3. 19) is known as the gas continuity" equation and is valid when 
applied to the system of gas species. Equation (3.20) is known as the particle 
continuity equation and is valid when applied to the system of particle species. 


The species continuity Eq.(3.6)may be written in a more convenient 
form by making use of the gas continuity equation and introducing the mass 
fraction of species i. In general, the mass fraction of species i may be 
defined as: 


X. 



(3.21) 


Substituting X . into the species continuity Eq. (3.6) 


(X. p ) +v # X- P q. - w. =0; i = l,NS 
8t i *rn' v i 'm n i 


expanding the first term 

dX . 8p 

i ' m — • 

P “ 5 r- +x. —sr~ +V • X P q. - w. = 0 i = 1 , NS 
r m 9t *1 ot i n i ’ 


(3.21a) 


and applying the vector identity 


V • X. p q. = V • (X.(p q.)) = VX- • P q. +X. (V» P q ) 

i m n i m n 1 m n i 
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results in 


m 


(ax. _ \ ia p 

TT + V v *d + ^TT 


+ V . p q 


i) - *i = 0 


i = 1 , NS (3.22) 


By limiting the system to only gas phase reactions (no reactions in the particulate 
phase) the species mass fraction X ^ may be redefined as: 


^ /p 

Assuming all gaseous species have velocity = q, Eq. (3.22) becomes 


(3.23) 


'("aT + V vx i^ + *i(!f +v« pq^- w. = 0 


i = 1, NS (3.24) 


Applying the gas continuity equation to the above expression the second term 
is set equal to zero and Eq. (3.24) becomes: 


/ax. \ 

’(sr + v vx i) - w i = 


= 0 i = 1, NS 


The final form of Eq. (3.24) is obtained by applying substantial derivative, 
Eq. (3.4) to yield 


DX. 

p = 0 i = 1,NS species continuity 

equation 


(3.25) 


Under the above assumption that there are no particle reactions, Eq. (2.20) 
for each particle species j becomes: 


+ V •p^q* = 0 j = 1,NP particle continuity 

equation 


(3.26) 
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For steady state flow processes, Eqs.*(3.13),(3.26)and (3.25) become: 

Equation (3.13) 

— & 

V ’P m < Jm = 0 global continuity (3.27) 

equation 


Equation (3.26) 

V . p j q j = 0 j = 1,NP particle continuity (3.28) 

equation, and 


Equation (3.25) 


# 

p q. • VX* - = 0 i = 1,NS species continuity (3.29) 

equation 


3.2.2 Momentum Equation 

The various forms of the global momentum equation for a chemically 
reacting gas-particle mixture are derived in this section. The forms of the 
equation derived, in order are: 


1. 

General Global Momentum Equation 

(3.42) 

2. 

General Global Momentum Equation 
with particle drag effects 

f 3.49) 

3. 

Global Momentum Equation for steady 
state, inviscid, no body force flow 

(3.50) 


The linear momentum of an element of mass m is a vector quantity 

defined as mq . The fundamental statement of Newton's law for an inertial 
m 

reference is given in terms of momentum as 

f ■ m (“mV 1 = 0 < 3 - 3 <» 

To derive the general global momentum equation the two terms of Eq. 
(3.30) will be analyzed separately. 
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The total force acting on the system may be represented as the sum 
of the surface force distributions (force distributions acting on the boundary 
of the system) and body force distributions (force distributions acting on the 
material inside the system) as follows: 


F 


NS _ 

= J a • n dS + J £ Pi f L dV 


(3.31) 


The surface force term may be converted into a volume integral by applying 
Gauss' theorem, Eq. (3.2). Therefore, 


F 



a dV + 



NS 


E p i f i dV 


(3.32) 


The surface stress tensor, a may be written in the form: 


a = -P6 + T 


(3.33) 


and 


V • a = - VP + V . t 


(3.34) 


Substituting Eq. (3 . 34) into Eq.(3.32) , the total force acting on the system 
may be written in the form 


F = 


NS 

f (- VP + V . r) dV + J £ p. t dV ; 


or, upon combining terms 

F = 


V i = l 


/(■ 


VP + V . r + 


NS \ 

E ^ ? i) dV 

i = l / 


(3.35) 
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The momentum term m m c l m ^n Eq.(3.30) may be rewritten in the 


form 


r NS 

"mV - J £ P . V dV ; 

V i = l 


from which the substantial derivative can then be written 


n n /* NS 

= mJ Z p i q i dV • 


NS 

z 

V i = l 

Applying Reynold's transport theorem Eq.(3.3)to Eq. (3.36) 

NS _ NS 


(3.36) 


5F - J Si Pi ?l dV + / V -E V V dV ; 

v i=l v i=l 


rearranging terms 


, NS NS , 

K |m mV = / If Z P . V + V ‘ Z P iV’i dV ; 

r, 1 ui i=i 


(3.37) 


and rewriting the results in terms of the gaseous and particle species results 
in: 


NP . „ 

lf>sq + V q^ 

l + V* (pqq+ y ?){dV 

II 

H- < 

V J =1 


or , 


Dt ,m mV 


V 


+ V • p q q + 


NP r . ni 

Zpt ^ A tV'pV <n| dV ( 3 . 38 ) 

j=l L J ’ 
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Substituting Eqs. (3.35) and (3.38) into Eq.(3.30), dropping the integral notation 
and rearranging the terms yields 


a - 


jj-(pq)+V.pqq + 


NP f a . . 

f E 9T 

j = l 


) +v 






NS 

-VP +V •T + '£p i P (3.39) 
i = l 


Expanding the first two terms of Eq. (3.39) 


■^(Pq) +v*Pq? = P^+q|f+pq‘Vq + qv*Pq; ( 3 . 39 a) 


at 1 H at 


and then applying the gas continuity Eq. (3.19) 


-^r + V * P q = 0 


8t 


to the results, yields 


0 —a _a — k 

8f(Pq) + V • P q q 


d a -» -*■ 

= P (-^- + q • vq) 


(3.40) 


The third term of Eq.(3.39) may be manipulated in a like manner using the 
particle continuity Eq.(3.26) to yield: 

Ellr (P j ?)+vp i ?t j l - + 

j=l l J j=l 


(3.41) 


Substituting Eqs. (3. 40) and (3.41) back into Eq.(3.39) results in: 

NP . _ NS 

-VP +V • r + £p. t ; 


+( 1 * V 


NP 

q) + E pJ (^" + ^ ,v ^) 


i = l 


or in terms of the substantial derivative, Eq. (3.4), 

_ NP . _ NS 

PdT« t iy 5T? - - VP + v.T* g Pi r. 
j=l i=l 


(3.42) 
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Equation (3.42) is the general global momentum equation. In this equation, 
the five terms are, respectively: (1) the non-stationary and convective rate 
of change of momentum per unit volume of the gaseous species; (2) the non- 
stationary and convective rate of change of momentum per unit volume of the 
particle species; (3) the net hydrostatic pressure force; (4) the viscous stress 
force acting on the surface of the unit volume; and (5) the body forces per unit 
volume. 

Since it is assumed that all forces are negligible, except that of the 
pressure of the gas and the particle drag, the only force exerted on the 
particle is a drag force caused by the relative motion between the gas 
and the particle. Figure 3-2 illustrates a spherical particle in a gas flow field. 



Fig. 3-2 - Particle' Momentum and Heat Transfer Model 

To determine the effects of particle drag on the global momentum 
equation, Newton's Second Law is applied to the drag forces acting on 
particle species j. 

Sj = 
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or rewritten in terms of particle density 


In general, the drag force on a spherical particle may be written in the form: 

d{ = C J D 7 T (r j ) 2 (|p A? )A?’|) 


where 


-x: — i ~x\ 

Aq J = q - q J . 


Equating the two expressions for particle drag 


ir(r^ p Aq^ | Aqj|; 


and rearranging terms yields: 


D -i 

n J 


3 (j c L\ 

8 W rV 


A? |A?| 


Equation (3.44) can be simplified by referencing to the Stokes flow regime 
(Reynolds number is less t. n 1) in the following manner. 


Define 


where 


Stokes 


; D \ = M ’ Nu = 2 

/Stokes R 


Therefore, 
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Substituting Eq.(3.45) into Eq.(3.44) and defining the parameter 


A i _ 9 

A 2 


V i } 

. .2 5 

m-* r-* 


the momentum balance becomes 


(3.46) 


mV - a! . 


(3.47) 


Expanding the substantial derivative by means of Eq.(3.4), Eq.(3.47) becomes 


q j + qj • V? = Aq^ = ^ • 


(3.48) 


Substituting Eq.(3.47) into the general global momentum equation results in: 

NP NS 


p|+^p ] A J A? = - VP + V • r + Pi f. 
j=l i=l 


(3.49) 


Equation (3.49) is the general global momentum equation with particle drag 
effects. 


For the case in which the flow may be described as steady state, inviscid 
and no body forces present Eq.(3.49) becomes: 

NP 

p q *Vq + 2 A ' J = * VP . (3.50) 

J = 1 

We note at this point that chemical reactions do not alter the forms of the 
global continuity equation or equations of motion. 

Equation (3.50) will be expanded now for later use in the derivation of 
the energy equation. For an arbitrary vector A the following identity exists 

A • VA = j V (A • A) - Ax f?x A) . (3.51) 
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Applying this ident.ty to the first term in Eq. (3.50) 

l -i —A — i — 1 

p q • Vq = yp v(q • q) - pq x(px q) ; 

substituting the results back into Eq. (3.50) 

N 

■| p V(q • q) - P q x(7 xq) + pj A J Aq^ = - VP ; 

j=l 


and rearranging terms yields 


NP 

\ V(q • q) = q x(V xq) - ± £ p j A j A? - 

j = l 


( .52) 


3.2.3 Energy Equation 

The various forms of the energy equation for a chemically reacting 
gas-particle mixture are derived in this section. The forms of ihe equation 
derived are: 


General Global Energy Equation Neglecting the 

Effects of 

(3 

.74) 

Radiation 





General Particle Energy Balance Equation 



(3. 

34) 

General Global Energy Equation with Radiation 

Effect s 

(3 

.85) 

Global Energy Equation with Radiation Effects 

for 

Plow 

(3 

. 105) 

Described as Steady State, Adiabatic. Inviscid, 

and 

No 



Body Forces Present. 






The most general form of the energy equation, according to the first law of 
thermodynamics or the law of conservation of energy can be wx itten as 


DE dQ dW 

Dt dt di 


(3.53) 
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The energy E in the control volume may be written in terms of the species 
present as 

-NS 

E = f ^e.dV; 


where 


V i = l 


q. q. 

e. = i. +-f +~Z. i = 1, NS 
i i 2 g c 1 


i. = internal energy of species i per unit mass 


(3.54) 


*i 


— 2 = kinetic energy of species i per unit mass 
Z i = potential energy of species i per unit mass. 


The substantial derivative can then be written 


DE 

Dt 


r> /*NS 

•&/Z> i*i dv 


NS 

E< 

V i = l 

and expanded by applying Eq.(3.3)to yield 


- NS NS 

V i=l V i=l 


(3.55) 


Recalling that the total force acting on the control volume according to Eq. 
(3.31) is 


F 


/ ■ •" dS + y £>. t dV ; 

S V i = l 


and applying the definition for the rate at which work is performed on or by 

, -a 

the control volume dW/dt = -q • F, 
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one obtains 


NS 


^ = - j*q . <T . n dS - J £ ?. . p. t dV 
S V i=l 

The surface integral may be converted into a volume integral by applying 
Gauss' theorem, Eq.(3.2). Therefore, 

HW /* = - r NS _ 

ST = - f Vto.q.)AV-f 

V V i=1 

Recalling the definition of the surface stress tensor o , Eq. (3.33) 


a - - P6 + T , 


and taking the dot product 


CT • n = - Pq + 7 • q , 

the first volume integral in the above equation may be written as 

J v • (- Pq + r . q) dV . 

V 

Substituting Eq.(3.56) back into the expression for dW/dt 

_ NS 

V / i = l 


combining integrals 

NS 

^ = .y[^.(-p?+?.5 + i;? i .p i r]dv, 

V i = l 

and expanding terms yields 

r _» =_ a> NS _ k 
= J [V • Pq - V • ( r . q ) - £ q. . p i t J dV . 


4W 

dt 


(3.56) 


(3.57) 
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The rate at which heat is transferred to or from the control volume may be 
written in the form 

33 =/S-?ds , 
s 

or in terms of a volume integral by applying Gauss' theorem, Eq.(2.2) 

f = / V.Q d V, (3.58) 

V 

where Q is the conduction heat flux vector. 

Substituting Eqs.(3.55),(3.57) and (3.58) into Eq.(3.53) and dropping the 
integral notation results in 

Ns NS ^ ^ _ _ NS _ 

5T ®i + V # S P i ®i *i " V * Q + V * ‘ V * (T **5) * E %. * P i *i = 0 * <3 * 

i=l i=l i=l 

Equation (3.59) is the general global energy equation neglecting the effects of 
radiation. This equation is based on a unit volume. The six terms are, re- 
spectively: (1) and (2) the total rate of increase of internal, kinetic and potential 
energy caused by local and convective changes; (3) the heat conducted to or 
from the control volume; (4) and (5) the work done on the fluid element due to 
surface forces; and (6) the work done on the element by the body forces. 

In deriving Eq.(3.59), all that has been neglected is the energy trans- 
ferred due to radiation, which if necessary, can be added to the equation as 
an extra term. Radiation effects will be determined later in this section. 

The microscopic quantum effects (interchange of energy and mass) are 
not considered in Eq.(3.59) as it has been implied by the classical first law 
of thermodynamics. The general global energy equation will now be expanded 
and written in terms of the gaseous and particle species present. 
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Recall that 2 __ 

q. q. 

e. = i. +- 5 -+—Z. i = 1 ,NS . (3.54) 

* 1 c 8 C 1 

If the potential energy term is combined with the body force term, Eq. (3.54) 
may be rewritten in the form 

2 

e. = i. +- 5 - i = 1 , NS (3.60) 

1 1 6 


Defining the specific enthalpy of species i as 

P. 

h. = i. +-*- i = 1, NS 
1 1 

where P. = partial pressure of species i 
and combining Eqs.( 3. 60) and (3. 61) yields 

P- q 2 

, i,i 

e. = h. - — + 

1 1 p. 2 

*i 

Next, define the total enthalpy per unit mass of species i as 

2 

% 

H. = h. +- 7 - . 

1 1 2 


(3.61) 


(3.62) 


Therefore, the expression for ei becomes 


e. 

1 


H. 

1 


P. 

1 



Substituting the expression for e. into Eq.(3.59) results in 



- V • Q +V • Pq - V • 


NS 

(T*q) VP/i = ° 


(3.62a) 
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The pressure P may be expressed as the sum of the partial pressures as 
follows NS 

P = . (3.63) 

i=l 

Substituting Eq.(3.63) into Eq. (3.62a) 

NS NS NS 


i=l 


i=l 


i=l 

NS 


+ V • Pq - V • (T • q) - ^ q i • pj fj = 0 

i=l 


and noting that 


NS 

V.£?i P i = v.pJ, 

i=l 


the above equation becomes 


NS 


NS 


NS 


^2>i H i - if + V •£ Pi ?i H i - v • b - v • ff • 5 - £ ^ • Pi h = 0 (3 - 64) 

i=l i=l i=l 


From the assumption that thermodynamic local equilibrium exists, the follow- 
ing expressions may be written in terms of the gaseous and particle species: 


and. 



NS 

E p i H i 

NP . 

= P H+ 

(3.65) 


i=l 

j=l 


NS 


NP . 


E 

p. q. H. = 
*i n i 

P q H + 2 P 3 q J n J 

(3.66) 

i=l 


j=i 



Substituting Eqs.(3.65) and (3.66) into Eq. (3.64) 
NP v NP 


( , \ NP . 

pH ♦ gpi - f + V . H i) - V. 3 - V M?-q) 


NS 

‘E V Oi 1 ! - °. 

i=l 
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and rearranging terms yields 


NP r i 

^(PH) + V*pqH+ £ O^H*) +V .p J ?H j = |f +V. (T.q) 

j=l 


NS ^ 

+ ^ q 4 . Pj f . + V • Q ; 
i=l 


where 


(3.67) 


H = h + q /2 (3.68) 

H J = h J + (q J ) /2 . (3.69) 


Expanding the first two terms of Eq. (3.67) 

P ff + H ^ + pq *V H + K V «p? . 
and rearranging terms yields 

p (f +?-™) + H (j£ + v. P s) . 

Recalling the gas continuity Eq. (3.19) 

+ V • P q = 0 , 

and applying the definition of the substantial derivative to the first term in 
the above expression 


9H - f-7 TT DH 
at + q . VH - , 


the first two terms of Eq. (3.67) reduce to 


|(P«) +V- pq H = p^ 


(3.70) 
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The third term of Eq.(3.67) may be manipulated in a like manner. 
Expanding 



VH^ + H^V • p 


j ?) 


and rearranging terms yields 



“5T + ^* 


VH J 




+ V 


P^ <? 



Recalling the particle continuity Eq. (3.26) 

+V • p j q* = 0 j = 1,NP , 

and applying the definition of the substantial derivative to the first term in 
the above expression 


dt q Dt ’ 


the third term of Eq.(3.67) reduces to 



(p j H j ) 


+ V • p^ qj 




j=l 


(3.71) 


Substituting Eqs.(3.70) and (3.71) into Eq. (3.67) results in 


P 


DH 

Dt 



PH* 

Dt 


~ + V 

at + 


(T • q) + £ ?i • Pi V • Q • 
i = l 


Recall that 


H j = h j + 


(3.72) 


(3.69) 
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Therefore, 


or, 


PH* 

Dt 


finally. 



(3.73) 


Substituting Eq.(3.73) into Eq. (3.72) results in 


P 


DH 

Dt 



+ 7 


Dt 


! ) 


NS 

ftp = -»■ -*■ -** 

gt + V • (T • q) + 2^ Qi* Pi f i +V * Q (3.74) 

i = l 


Equation (3.74) is the general global energy equation written in terms of the 
gaseous and particle species. 


To determine the effect of particle/gas conduction and radiation, the 
particle energy balance for the j** 1 particle is written as follows: 


|4jt (r j ) 2 j (T j - T) -o4» (r j ) 2 (T j ) 4 - 0t j T 4 J ; (3.75) 


where 


. k 

= particle heat transfer 
2r^ film coefficient. 


(3.76) 


Defining the Nusselt number parameter 

= 

where 


N J 

u 


(Nu) 


(Nu) 


Stokes 


Stokes 
= 2 ; 


(3.77) 


(3.78) 


3-25 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC TR D867400-I 


the particle heat transfer film coefficient may be rewritten as 

J 


j - *2: 


K J = 


• • 




(3.79) 


Furthermore the heat quantity of the particle is defined as 

. 3 


Differentiating Eq. (3.80) 


rearranging terms 


Q j = (rJ) h J. 


20? = / r i) 3 Shi 

Dt 3 * ' ’ Dt ’ 


(3.80) 


Dh J 

Dt 


DQ J 


3 Dt * 

4jr m J (r J ) 


and substituting the above results, with the expression for back into Eq.(3.75), 
the particie energy balance equation becomes 

m j^ (r J, 3 ^ . . !^i | 4 „ ( r l) 2 j ( T j . T) - ,ri> 2 Ut^T 4 ). 

Solving for Dh^/Dt in the above equation 

^ = - k -4 ^ <Tj - T) - -TT k j(Ti)4 - “ j T 1 • 

(rJ) m J m J r J l 1 


and recalling Eq. (3.46) 


A j _ 9 -Jdi. 

A - z . .2 

j i 

m J r J 


(3.46) 


the equation for Dh^/Dt may be written as 


^ ■ - % sty < t3 - t > - -fr [*W ■ *i 

V f J m J r J 1 1 


(3.82) 
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Equation (3.82) may be simplified by letting 

C j = 2 _!e. = (3.8: 

r Pr i/O 

Therefore, 

= - \ pj C * (T-> - T) (T^) - T 4 ] . (3.84 

m J r J l J 

This is the general particle energy balance equation. The temperature of the 
particle depends upon its state, where the state may be a liquid, a liquid in the 
process of solidifying, or a solid. The temperature is uniquely related to the 
particle enthalpy by the equation of state T^ = f(h^), tabulated. 


Substituting Eqs. (3.47)and (3.84) into the general global energy Eq. (3.74) 


yields 


NP I 

p TfiT + 2 - f p j A j C j (T j - T) - Y j [£ j (T j ) - ot j T 4 J + q* • p j Aqj 

j=l m F 

NS 

= H+V .(T.q) • Pi t + V • Q . (3. 


This is the general global energy equation with radiation effects. 

Expanding the first term of the above equation 

_ DH „ /9H . -* 
p Dt ~ p (at + q * VH ) * 

and applying the definition for total specific enthalpy 


H = h +3- ; 
2 


the first term becomes 


DH AH -» I a a 

p^kt =p^r + pq ,vh+ p < i , 7 y (q*q) • 
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Recalling the momentum balance relationship 

N p 

\ V(q • q) = 5 x (vxq) - ^ Aq j - ~ , 

j=l 

and substituting it into the above equation yields 

NP 

P 15t = P Itf + P *** ( Vh ' + P ?• [7* xqj] - if. jjj p^ Aq^ . 

j=l 

This expression may be further simplified by noting that for an arbitrary 

«ik 

vector A the following identity exists! 

A'(Ax^xA)] = 0. 


Applying the above identity to Eq. (3.86) 
DH JH, - 


• Pf ♦Pq-* h -^>-q.£p'A>*? . 

j=l 

and substituting the result back into Eq.(3.85) yields 

NP NP 

p w + p ^ • <vh ■ ~ p~> * y! pj ^ 

j=i j=i 


- | p j A j C j (T-! - T) 


-i£_. p i gi (T J) . a j T 4 ] + J . p j A j Aq j 


an j, 

= gj +V. ( T .q) 


NS 


+V q| • + v • Q . 

1^1 

For the case in which the flow may be described as steady state, 
adiabatic, inviscid and no body forces present the above equation becomes 

NP . NP i 

pq . (Vh - ~) - q . p J A J A? + ^ j- § p^ (T^ - T) 

j=l j=l I 

r 2 -? p j U j (T j ) - aJ T 4 j + q j • p j A j A? 

m^ r^ L * 


= 0 
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This equation can be reduced t a simpler form by letting 


B j = q . Aqj - P . + \ C 3 (T j - T) + 


3 q 

• ♦ • 
m-* r* 



(3.86c) 


Therefore, the global energy equation becomes 

NP 

Pq ‘(Vh - —■) - T]p J A J B j = 0 (3.87) 

P j=l 

To expand the term in parenthesis in the above equation, it is necessary to 
examine the fundamental thermodynamic relationship for a system in which 
chemical reaci.ons are occurring. 


3.3 GASEOUS THERMODYNAMIC RELATIONS 

Assumption 4 given in Section 3.1 states that the gas obeys the perfect 
gas law and is in chemical equilibrium, nonequilibrium or is chemically 
frozen. A chemically frozen gas is one in which the gas stops reacting at 
a given species concentration so that the molecular weight along a stream- 
line remains at a fixed value and the ratio of specific heats is a function of 
temperature only. However, for the chemically reacting case the local 
chemical species concentrations will change in accordance with type of 
chemistry assumption considered for a respective analysis. 


In the high-temperature low-velocity regions of the flow field, an equi- 
librium chemistry calculation for the gas phase is a good assumption. Local 
residence times of the flow are sufficiently long for all chemical reactions to 
proceed to completion. However, as the flow accelerates through the nozzle 
and exhaust plume the local flow residence time becomes less than that re- 
quired for the chemical reactions to reach completion. Significant deviations 
from the chemical equilibrium conditions occur and ultimately the flowfield 
chemistry usually approaches a frozen condition. This calculation is treated 
best by including the kinetics in the gasdynamic calculation, thereby avoiding 
the problem of choosing either an equilibrium or frozen chemical analysis. 
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The development of the gas thermodynamic relations and the corre- 
sponding contribution to the gasdynamic relations will consider the chemical 
kinetics. Appendix C addresses the equilibrium chemistry analysis for gas- 
particle flows. Appendix E discusses the nonequilibrium chemistry analysis. 
Finally, a summary of the applicable equations for chemical nonequilibrium 
and chemical equilibrium flow is presented in Section 4.2.3. 

Consider only the gas phase portion of the flow system in which there 
are NS-NP different species present. If the mass fraction of each species i 
is constant, the specific enthalpy, h, of the gas depends only on specific entropy, 
S, and pressure, P. However, for a variable composition 


h 


b (S, P, Xj» X^> • • • ’Xj^s-NP^ 


and thus the total differential of h is 


Hh _ 9h 

dh " as 


P.X ; 


oo . 8h 

dS + ap 


NS-NP 
dP + 


S.* 




dX 


In this expression the subscript X ^ implies that the mass fractions of all 
species are constant during the variation in question. On the other hand, 
the last term in the equation is a sum of partials in each of which the S and 
P are constant, together with all but one of the mass fractions. 


For constant mass fractions the total differential of h may be written 
as 

dh = T dS + ^ dP . 

r 


Therefore, 


9h 

as; 


P.*; 


T and 


ah 

ap 


s.x, 


i 

p 


By defining the chemical potential Mi as 


Mi = 


ah 

ax, 


= h. 


i’s, P,Xj 


specific enthalpy of species i 
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the fundamental thermodynamic relationship for a system in which chemical 
reactions are occurring may be written in the form 

NS-NP 

T <**1 (3.88) 

i = l 

NS-NP 

y: m i vx L 0.89) 

i = l 

Since the pressure is a function of the gas state variables (p, S) 

P = P(p, S) , (3.90) 


dh = T dS + ^ + 

r 

The above equation may be rewritten as 

VP 

Vh - = TVS + 


the expression for VP may be written 

vp • (§) s vp+ (f|, vs 

Noting that 


and. 


2 _ /0P\ 

a "U) s 

p Bridgeman's Equation 


(3.91) 


(3.92) 


7 P P 

where a, c^ and R are local thermodynamic equilibrium properties, VP 
becomes 

VP = a 2 VP + — 5 VS . 

P ' 

Rearranging terms in the above expression 


VS 




(VP - a Vp) . 


(3.93) 
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and utilizing the equation of state for an ideal gas 


P = pRT , 


Eq. (3.93) becomes: 


VS 



(VP - a 2 


Vp) . 


Multiplying through by T 


Tvs = (~pir~) ^ p - a 2 vp) . 

and rearranging terms yields 

ws ■(tF- , )(?-^ b ) • 


Substituting Eq.(3.97) into Eq.(3.89) results in 


Vh 


VP 

P 


!e . A(VP _ a 2 y fi 

R 7\ P 0 


v NS-NP 

) + E 

i=l 


(3.94) 


(3.35) 


(3.96) 


(3.97) 


(3.98) 


Since particle si cies are not considered in the chemical reactions, the 
summation term upper limit may be defined as 

NG = NS - NP (3.99) 

where NG = number of gas species present. 


Equation (3.98) then Becomes 


7 NG 

i = l 


(3.100) 
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Substituting Eq.(3.100) into the global energy £q.(3.87) 


pq 


[(f - ') (ir - —^) + E "i v *i] - E p j B j = » 


and expanding the dot product results in 



^ NG ^ NP . . . 

(q . VP - q • a VP) + P £ q • VX. - ^ p J A J B J 

i=l j=l 


0 (3.101) 


Equation (3.101) is the expanded form of the global energy equation for flow 
described as steady state, adiabatic, inviscid and no body forces present. 


Dividing through by (c ^/R - 1) yields 

~ . NG DX. . NP . 

q.VP-a q » VP + c - /r ~ -W - c /R - I L. P* aJ fiJ 

P i=l P j=l 

Recalling the species continuity equation 

DX. 

P = w. i = 1, NS , 


= 0. (3.102) 


(3.25) 


and letting, 


E i B J 

B 1 " r7RTT 


(3.103) 


and 


^i " rjR 

p x i=i 


NG 

br2 | ' 1 w i ■ 


(3.104) 


the final form of Eq.(3.102) becomes 


NP 


q.VP-a q. Vp + p J B j = 0 

Pi 


(3.105) 
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3.4 SUMMARY OF THE GOVERNING EQUATIONS FOR STEADY, ADIABATIC, 
NONEQUILIBRIUM FLOWS OF REACTING GAS-PARTICLE MIXTURES 
WITHOUT TRANSPORT OR BODY FORCE EFFECTS 


The system of basic governing equations for any reacting flow field have 
now been derived. When the reactions between components of the gas mixture 
are known and the boundary conditions adequately specified, one should be able 
to solve the nonequilibrium flow system. 


The pertinent equations derived in this section are summarized as 
follows: 


Continuity equation for steady state flow 


V • p q = 0 

gas continuity equation 

(3.19) 

V*P j ? = 0 j = 1,NP 

particle continuity equation 

(3.26) 

p q. VX. - w. = 0 i = l.NS 

r ^ X 1 

species continuity equation 

(3.25) 


Momentum equation for steady state, inviscid flow with no body forces 
NP 

pq.yq+T" p J A J Aq J + VP = 0 global (gas + particle) momentum (3.50) 

• • - 5 Ami ^ 1 1 


j=l 


equation 


where 


A j _ 2 f-JdL 
A ~ 2 . .2 
LmJ rJ 


(3.46) 


p j 5.v? - p^ A^Vq^ =0 j = 1,NP particle momentum equation (3.47) 


Energy equation for steady state, inviscid flow, with no body forces and no heat 
loss from the gas-particle system 

NP 

q . VP - a 2 q • VP - T* P^ A^ B^+ 4*, = 0 global energy equation (3.105) 

FI 
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where 


A % 

B 1 = r /A 1 ^ + f cj < Tj - T > + -j ^ j [^ j (T j )’ -o j T 4 J <3 


A? = q -? 


and 


d = 


k G j 


qj . Vh j + | A j C j (T j - T) + -4 2 -t I (T j ) 4 - a j T 4 ] = 

r l ' 


.103) 


(3.43) 


(3.83) 


j = 1, NP (3. 

particle energy - 
ba lance equation 
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4. THE METHOD OF CHARACTERISTICS SOLUTION 

The set of equations summarized in Section 3.4 are applicable for flow 
calculations in both the subsonic and supersonic regimes. Examination of the 
equations reveals that a closed form solution is not possible, thus necessitating 
a numerical solution to obtain the flowfield properties. The choice of the 
numerical solution is governed by the application and the flow regime of 
interest. 

The current study addresses the nozzle-plume flowfield where the flow 
is everywhere supersonic except for those regions in the immediate vicinity 
downstream of a normal shock. Treatment of these regions are considered 
beyond the scope of this study, thus reducing the problem to one of analyzing 
a supersonic flowfield. For supersonic flow applications the method-of- 
characteristics provides a reliable method for numerically solving the 
set of governing equations. The method is characterized by rapid convergence 
of the numerical solution and has been shown to be unconditionally stable. This 
method is the one chosen for the current study. 

A number of investigators (Refs. 8,9 and 10) have developed solutions 
employing the method-of-characteristics. In "characteristic form" the 
relations are transformed to a set of differential equations which apply along 
the characteristic directions. Prozan (Ref. 16) developed a nozzle-plume: 
solution which includes gas thermochemistry considerations. Kliegel (Ref. 9) 
and Hoifman (Ref. 8) extended the method-of-characteristics to include the 
treatment of gas -particle flows. The present study followed the work of 
Prozan in the treatment of the gas equilibrium thermochemistry considera- 
tions anc that of Thoenes, etc., (Ref. 39) for chemical kinetics. These con- 
siderations have been combined with the coupled gas -particle solution as 
formulated by Kliegel, 
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The original formulation by Prozan utilized the enthalpy-entropy- 
velocitv form of the onmn»Mhii»ty equations, i.c., independent variables «*re 
the velocity, total enthalpy and entropy. Prozan showed that the thermo- 
chemistry calculations could be uncoupled from the gasdynamic solution. 
Consequently the thermodynamics were calculated in terms of the independent 
variables, V, H, S and retrieved via interpolation from tables as needed by the 
characteristic solution. The routines became an integral part of the stream- 
line normal code developed by Ruo (Ref. 17) from which the RAMP code evolved. 
Consequently the V, H, S form is used in »>.e present analysis for chemical 
equilibrium calculations. However, when kinetics or transition between the 
continuum or non-continuum flow is to be treated, use of the pressure-density- 
velocity form of compatibility equations is more appropriate. Therefore 
development of both forms is given in the following sections. 

4.1 DEVELOPMENT OF THE CHARACTERISTIC CURVES 

The governing flow equations summarized in Section 3.4 may be written 
in the following expanded two-uimensional (or axisymmetric) form: 

Continuity equation for steady state, particle species not considered in the 
chemical reactions 


(3.19) pu +pv + up + v p -f s 0 

x y *x *y y 

(3.26) p^u^ + p V + uV + vW + v = 0 j - l.NP 

’ x y x y y J 

(3.25) pu (X.) +p v (Xj) - w. = 0 i = 1, NG 

x y 

Momentum equation for steady state, inviscid flow, with no body forces 

NP 


(3.50) 


puu„ + flvu + P 
x r y x 


+ A^ (u - u^) 

j=l 


0 


(4.1) 


(4.2) 


(4.3) 


(4.4) 
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nnv +Ovv +P + V' (v - v^) = 

* X y y 


(3.47) 


ir* + p^ - p* (u - u^) = 0 j = 1 , N P 
x y 


(4.6) 


p* + p 3 v** - p 3 (v - v^) = 0 j s 1 , N P 

x y 


(4.7) 


Energy equation for steady state, inviscid flow, with no body forces and no 
heat loss from the gas -particle system 


(2 105) u P x + v P y - a u P x - a v P y " ^ p J A J B| + = 0 (4.8) 

j-1 


( 3 . 84 ^ 


pVh j +pVh j + | p j A j C j (T j - T) + (T j ) 4 - T 4 ] =0 j= 1,NP (4.9) 

x y m 3 r 3 L J 


The corresponding coordinate system used in this analysis is presented 
in Fig . 4-1 . 



Fig. 4-1 - Coordinate System Used in the Development of the Characteristic 
and Compatibility Equations 
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To minimize the amount of "bookkeeping" that would be required in the 
development of the characteristic equations the governing flow Eqs. (4.1) 
through (4.9)will be written in the following general form 

8 * d<t> 

L m = a mn + u mn W + c m for n = and m = 1,N (4.10) 

where: 

N = 4 + 4NP +NG 

* ♦ * t 

♦ n represents the dependent variables u, v, P,p,u , v^,p^, h^ andX'. 

a and b are coefficients. 

mn mn 

If we let 

=u 

*2 =v 
* 3 =P 

* 4 =P 

*4J+1 »> J.l.NP 

*4j+2 = TJ . J=llNP 

*4j+3 =pJ . j = 1 * NP 
^ 4j+ 4 =hj j = 1, NP 

*4 + 4NP +i = X i { ~ 1,NG; 
equation (4.10) can be rewritten as 



L = a u + b , u + a , v +b ,v +a - P +b , P 
m ml x ml y m2 x m2 y m3 x m3 y 

+ a a P 4b . p +a , u* + b _u^+a ,v^+b , v ^ 
m4 m4 H y m5 x m5 y m6 x m6 y 


+ a m7 p x + b m7 P y * a m8 h x 


+ a a h* + b 0 h^+.. 
m8 x m8 y 


+ a m, 4j+l u i + b m , 4j+l u y + * * * + a m, 4j+4 h x + b m, 4j+4 y + * * 


+ a m, 4 +4NP +i ^i^x b m, 4 +4 NP +i ^i^ + c for m = 1, 


y m 

y m=l 
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We then have N independent, linear nonhomogeneous equations written in N 
unknown derivatives. The N linear equations for may be combined to 
form a single differential operator by employing arbitrary multipliers and 
summing. Thus, 


L 


N 

= jr a L 

m m 

m=I 


0 


(4.13) 


where = arbitrary multipliers. 


Assuming that the following relation holds. 


^ = b cr for n = 1,N and m = 1 , N ; 


3. (J j - w v 

mn m dx mn m 


(4.14) 


Eq. (4.13) can be rewritten in the form of an ordinary differential equation. 
Substitution of Eqs.(4.10) and (4.14) into Eq.(4.13) yields 


o.* 


N 8<l> . 8$ 

L - V' (a a -jr~ + cr a -rjy -k— + c o ) = 0 
/ f m mn ox m mn dx o y mm 

m=l 

N JM 8 <t> 

— 5 — dx + <7 a ■ a — _ 
m mn ox m mn oy ‘ m ;n 


L = V (<y a -jp-S dx + a a - a ” dy + c a dx) = 0 

' tvi mrt )lv wi rv»rt H ir ' m • n 


m = l 


and finally, 


N 


L = ^ (a a d0 + c o' dx) = 0 n = 1,N . (4.15) 

u m mn n mm 

m=l 

Equation (4.15) is the generalized compatibility equation and is valid if and 
only if Eq. (4. 14) is true. 

Equation (4. 14) may be rewritten in the form 


O' (a ^ - b ) = 0 n = 1 , N ind m = 1 , N . 
m mn dx mn 


(4.16) 
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If Eq. (4.16)has a solution other than trivial, i.e., all a m = 0, then the 

determinant of the coefficients of a must be zero. Therefore: 

m 


D 


a b 

mn dx mn 


0 n = 1,N and m = 1,N 


(4.17) 


The equation D = 0 is called the characteristic equation for the system of 
equations (4.16). On expanding the determinant, it is seen that Eq. (4.17) 
is an algebraic equation of degree N and thus has N real roots. These 
roots are called the characteristic roots. 


To analyze the determinant D, begin by substituting Eqs.(4.1) t (4.4), 
(4.5)and (4.8) into Eq. (4. 12), and putting the results into a matrix of the 
form: 


a ll’ b ll 

a 12' b 12 

a 1 3 * b 1 3 

a 14' b 14 

a 21’ b 21 

a 22’ b 22 

a 23 ,b 23 

a 24’ b 24 

a 31’ b 31 

a 32’ b 32 

a 33 ,b 33 

a 34* b 34 

a 4 1 ’ b 4 1 

a 42’ b 42 

a 43’ b 43 

a 44' b 44 


This yields 


P, o 

o.p 

0,0 

u, V 

pu.pv 

0,0 

l.o 

0,0 

0,0 

pu.pv 

0, 1 

0,0 

0,0 

0,0 

u, V 

2 

-a u 


and 



NP . 

c 2 = ^ P Al (U ' " j) 
j=l 


(4.18) 


(4.19) 


(4.20) 
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NP 

C, = E P j A j 

j=l 


(v - v J ) 


NP 

C 4 = " E P 3 B 1 

j=l 


(4.20) 

Cont'd 


Next, substituting Eqs.(4.2), (4.6), (4.7) and (4.9) into Eq. (4.12) with j = 1 
(for the first particle) and putting the results into a matrix of the form of (3.18) 
yields 


where 


p \ 0 

pV.pV 

0 . p ‘ 

u 1 , V 1 

0,0 


0,0 

0,0 

0,0 


0,0 

n 1 LI 1 
P u ,p V 

0,0 

0,0 

(4.21) 

0,0 

0,0 

0,0 

n l 1 1 1 

P u ,P V 


r 

_ 6p*v l 


1 


5 

y 

u 1 ) 



C 6 

= -pWfu- 


(4.22) 

C 7 

= Va 1 ( v - 

V 1 , 

J 



C 8 * 


= IaWcT^-T) + -42, [^(tV-«^T 4 1 

r J l J 


(4.23) 


Repeating for each particle, j, and generalizing yields 


p\, o 

0 . p j 

J,v> 

0,0 



P J U J , p J V J 

0,0 

0,0 

0,0 



0,0 

pU.pV 

0,0 

0,0 

j = 1 , NP 

(4.24) 

0,0 

0,0 

0,0 

P J U J , P J V J 
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and 


4j+l 

C 4j+2 

°4j+3 

C 4j+4 


. LQ. 




- (u - 4) 

- p^ (v - v^) 

* 

i 


p^ 4 


j = 1, NP 


Substituting Eq. (4. 3) into Eq. (4.12) 


a m,4+4NP+i = 0 

if 

m ^n 

= pu 

if 

m = n 

b m, 4+4NP-K = 0 

if 

m \n 

= pv 

if 

m = n 

C 4+4NP+i = ' W i 

i = 

1,NG; 


and rewriting in the matrix form of (4.18) yields 


n = 4 +4N P+i 


pu.pv 

0, 0 


0 , 0 


0,0 0,0 


pu,pv 

0 ,0 


0 , 0 
• pu,pv 


If all of the coefficients of Eq. (4.12) were put into a single matrix of the 


a U ,b ll a l2 ,b l. a IN ’ b IN 

a 21* b 21 

X = 1 

a Nl ,b Nl a NN ,b NN 
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(4.25) 

(4.26) 

(4.27) 

(4.28) 
form: 

(4.29) 
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where 

N = 4 + 4NP +NG ; 
the result would be the diagonal matrix 


\ = 


*1 

0 

0 . 

. » . 0 • * • • 

0 

0 

0 

X 2 

0 . 

• i i 0 # * # » 

0 

0 

0 

0 

*1- 

... 0 . . . . 

0 

0 

0 

0 

0 . 

• . . o . . . . 

0 

0 

0 

0 

0 . 

J 

• • • A 2 • • • • 

0 

0 

0 

0 

0 . 

• 0 • • • • 

0 

0 

0 

0 

0 . 

... 0 . . . . 

,NP 

A 2 

0 

0 

9 

0 • 

■ 0.... 

0 

\ 3 


(4.30) 


Examining matrix (4.30), from which the coefficients a and b for 

mu mn 

determinant D will come, it can be seen that D will have the same form as 
matrix (4.30). That is: 


where 


D = 


D 1 

0 

0 . . . 

. 0 . 

... 0 

0 

0 

D 2 

0 . . . 

. 0 . 

... 0 

0 

0 

0 

-> 

D r * • 

. 0 . 

. . . 0 

0 

0 

0 

o . . . 

. 0 . 

... 0 

0 

0 

0 

0 . . . 

■ D 2‘ 

... 0 

0 

0 

0 

0 . . . 

. 0 . 

... 0 

0 





_NP 


0 

0 

0 . . . 

. 0 . 

D 

• * ' u 2 

0 

0 

0 

0 

0 . 

... 0 

D 3 


ji>j| obtains its elements from 
from etc. 
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For purposes of simplification, let: 


S = u{j£-v 

si . u i 4c . vJ 

dx 

y ' = *Y- 
Y dx 


(4.32) 


Substituting the coefficients from matrix (4.19) into matrix (4.17) and using 
Eq. (4.32)results in 


py' 

ps 

0 

0 


-P 

0 

ps 

0 


0 

y' 

-l 

s 


s 

0 


-a^S 


(4.33) 


Matrix (4.33) may be rewritten in the form 

N ■ ■ 


where 

and 


j = l 


d. . = elements of D. 
iJ 1 

N. . = the minor of d. . . 
ij U 


(4.34) 


Utilizing Eq. (4.34) and letting i = 3, Eq. (4.33) becomes 

' # 


1 


= ^(-D 4 n 


,0 


31 + d 32 N 32 + d 33 (-1)6 N 33 + >34 ( " 1)? N 34 * 


Expanding terms 



py' 

0 

s 


py’ -P s 

= -pS 

ps 

y' 

0 

+ (-l) 

p3 0 0 


0 

s 

-a 2 S 


0 0 -a 2 S 
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performing the necessary matrix algebra 


|dJ = -PS [-pa 2 S(y') 2 +pS 3 ]-l (-p 2 a 2 s 2 ) ; 

and combining terms yields 


1 


= p 2 S' |a 2 (y') 2 -S 2 + a 2 ] . 


Substituting the expressions for S and y' from Eqs. (4.32) 

i2 


D i| = p2 (" S - v ) 2 [ a 2 (S) - ( u ^ - v ) 2 + aZ ] ■ 


expanding 

Id 


= P 2 (" S' ■ 'V | a2 (S) ' " 2 (S) +2»v^-v 2 +a 2 j; 


and recombining terms yields the final form of |D^ 

'ii = p z ( u & ■ v ) 2 [ (a2 • u2) (S) + 2 uv S +(a2 


v 2 )|. 


(4.35) 


The game procedure is applied to matrix (4.24). Substituting the coefficiei 
from matrix (4.24) into matrix (4.!7)and using Eqs. (4. 32) results in 


lD^I = 


pV 


J ci 


P y o o 

0 P& 0 

0 0 0 


0 

0 

0 

P y 


j - 1,1MP . 


(4.36) 


Utilising Eq.(4.34)and letting i = 2, Eq.(4.36) becomes 

Wt\ - d 21 <- 1,3n 21 '•‘» 4n 22 <- 1)5n 23 + 4* I' 1 ' 6 * 
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Expanding terms, performing the necessary matrix algebra 

-p j S* 0 

p j 0 0 


mi = p j s j (-D 


and combining terms yield 


0 pj s J 


.2 .3 

p j s J [-(pj) (S j ) ] , 


ml = ot > j ) 3 (s j ) 4 


Substituting the expression for from Eq.(4.32) yields the final form of |D^ 

3 - - 


HI - ^ • 


2* 

(4.37) 


Finally, the same procedure is applied to matrix (4.28). Substituting the 
coefficients from matrix (4.28) into matrix (4. 17) and using Eq. (4.32) results 
in 



(PS) 


NG 


Substituting the expression for S from 7! j. (4.32) yields 



(4.38) 


Since the -aatrices |Dj| through | | are the diagonal elements of matrix |D[, 
the characteristic equation for the system of Eq.(4.16) can be written as 


Dj * 


D. 


* DJ 


NP 

D z * 


= 0 


(4.3C, 
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Substituting Eqs.(4.35), (4.37) and(4.38)into Eq. (4.39) 

v2 


|D| = P 2 (» f* - v) 2 Jfa 2 - u 2 ) + 2 uv ^ + (a 2 - v 2 )J p NG (a v) 

(pi) 3 (ai^-vJ) 4 . 0 . 


and combining termr • the final fom: of the characteristic equation becomes: 

v 2 . > >1.3/.. .v4 


,NG 


l D l"( u to- v ) | (a2 - u *> (S) + 2 uv * + (a ‘ 


] <P j »VsM’ -o 


(4. 


Realizing that p 1 ^ o, the characteristic roots of Eq.(4.40) are determined 
by setting each of the three remaining terms to zero and solving for dy/dx. 


Setting the first term to zero 



and solving for dy/dx results in 


djr 

dx 


— = tan9 . 
u 


Setting the fourth term to zero 



j = 1 , N P 


and solving for dy/dx results in 



dx 


X1 J 


tanQ^ 1 


j = It NP . 


(4.41) 


(4.42) 


Referring to Fig. 4-1, Eqs. (4.41) and (4.42) shows that each of the character- 
is'.aCs is inclined at an angle tangent to the gas velocity vector and particle 
velocity vectors re-^ectively. Thus we see that the characteristics are 
identical with the b j streamlines and particle streamlines of the flow. 
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The final characteristic root is obtained by setting the second term in Eq. <4.40) 


to zero 


(a 2 - u 2 ) + 2uv ^ + (a 2 - v 2 ) = 0 , 


and solving for dy/dx in the following manner: 


Applying the quadratic formula 


djr _ 

dx 


- 2 uv + ^4 u 2 v 2 - 4 (a 2 - u^) (a 2 - v 2 ) 
2 (a 2 - u 2 ) 


expanding the terms under the square root 


±L = 

dx 


-ZuvjZ Ku 2 v 2 - a 4 ± a 2 v 2 + a 2 u 2 - u 2 v 2 
2 (a 2 - u 2 ) 


dividing through by 2 and combining terms under the square root 

, J 2 7T~~T T 

dy - UV + ^|a (u + V ) - a 

dx = 2 2 

a - u 


and finally moving a outside the square root yields 


dx 


- uv + a 2 ^(u 2 + v 2 )/a 2 - 1 

2 - 2 
a u 


(4.43) 


The first term in the square root may be simplified by utilizing the definition 
of the gas Mach number 

M = q/a 


, Z 2 2 

m 2 = s_ , a+v. 


2 


(4.44) 


Substituting Eq. (4.44) into Eq. (4.43) yields 


= 

dx 


uv + a 2 ^M 2 - 1 

2 2 
a - u 


(4.45) 
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Referring to Fig. 4-1, the following relations may be written: 


and 


a = 

q cos0 

(4.46) 

v = 

q sin0 

(4.47) 

a = 

sin (l/M) = sin ”^(a/q) 

(4.48) 

a = 

q sina 

(4.49) 

cota 

= Jm z - 1 

(4.50) 


Substituting Eqs. (4.46) through (4.50) into Eq. (4.45) 


_ 

dx 


2 2 2 
- q cos0 sin9 q sin a cot a 

2 " . 2 27 

q (sin a - cos 0) 


dividing through by q 


_ 

dx ' 


-cosG sin0 + sina cosa 

2 2 
sin a - cos 0 


and multiplying the numerator and denominator by - 1 yields 


dy cos9 sin9 + sina cosa 

dx = 2. . 2 

cos 8 - sin a 

To simplify Eq. (4.51) consider the following trigonometric identities: 

sin(9 + a) = sin© cosa + cos0 sina , 
cos (9 + a) = co80 cosa + sin8 sina , 
cos^9 - sin^a = cos(8 + a) cos(0-a) . 


(4.51) 


(4.52) 

(4.53) 

(4.54) 


Multiplying Eq.(4.52) by (4.53) 


— 2—2 
sin(9 + a) cos (0 +a) = cos a sin0 cos0 + sin 9 sina cosa 

— 2 2 
+ cos 0 sina cosa .+ sin a sin0 cos0 


and combining terms yield 


— 2 2 
sin(9 + a) cos{9 + a) = (sin a + cos a) sin0 cos9 

— 2 2 
+ (sin 0 + cos 0) sina cosa ; 
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or 


sin(0 +a) cos(0 +a) = cos© sin0 + sina cosa . (4.55) 


By substituting Eqs.(4.54) and (4.55) into Eq. (4.51) the final form of dy/dx is 
obtained 


= 

dx 


sin (9 + o) cos (9 + a) 
cos(0 + a) cos (9 - a) 


sin (9 -f a) 
cos (0 ^ a) * 


or 


dx 


tan(0 + a) . 


(4.56) 


Referring again to Fig. 4-1 , Eq. (4.56) shows that each of the two char- 
acteristics is inclined at the Mach angle to the gas velocity vector. Thus we 
see that the characteristics are identical with the gas Mach lines of the flow. 


In summary, the characteristics of the flow have been found to be the 
gas streamlines, the particle streamlines (one for each particle species) and 
the gas Mach lines. 


4.2 DEVELOPMENT OF THE COMPATIBILITY EQUATIONS FOR GAS- 
PARTICLE FLOWS 


The relationship assumed in Eq. (4.14) 

^ = b <j for n = 1,N and m = 1,N 


a a , 
n.i m dx 


CT 

mn m 


is valid along each of the characteristics. Therefore, we can now substitute 
the governing flow equations into Eq.(4.14) and knowing the conditions that 
exist along each characteristic, solve for the multipliers, <y . The corre- 
sponding compatibility equations that are valid along each characteristic are 
then derived by substituting the proper multipliers into the general compati- 
bility equation. The compatibility equations will be first derived in the 
pressure-density-velocity form. This form of the equations is used by the 
RAMP computer program when the finite rate chemistry option is utilized. 

In this form, the local temperature may be solved for directly eliminating 
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the need for a time consuming iterative solution. The compatibility equa- 
tions will then be rederived in the enthalpy-entropy -velocity form. Tb’S 
form of the equations is used by the RAMP computer program when the 
equilibrium and/or frozen chemistry option is utilized. In this form, the 
thermochemical data generated by the TRAN 72 computer program may be 
used. 


4.2.1 Pressure-Density-Velocity Form of the Compatibility Equations 

Substituting the governing flow, Eqs. (4.1) through (4.9) into Eq. (4.14) 
results in 

for n = 1: 

ff i p 5^ + ff 2 < pu *F- pv > = 0 

or 

a iS +a z (u i£- v) = 0; 

for n = 2: 

- <*1 P + P (u - V) a 3 = 0 

or 

°1 * (a to - v) “5 ! 

for n = 3: 

(T 2^- <T 3 + (U S- V)<T 4 = 0; (4.59) 

for n = 4: 

( U ■ v ) a l ■ &2 ix ' V ) a 4 " 0 » (4.60) 

for n = 4j+l: 

pj to'4j+l +p) ,aJ to * v) > 0 4j+2 = 0 J ' l.NPi 

or 

to a 4j+l + ( ' jJ to ' vJ) a 4j+2 = 0 J ■ >' NP ! 0.41) 
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for n = 4j+2: 


P J ^4j+l + pJ (uJ to ’ vJ) °4j+3 = 0 j = ! ’ NP 


or 


°4j+l * <" J I* - ^ d 4j+3 j = l.NP ; 


(4.62) 


for n = 4j+3: 


for n = 4j+4: 


or 


fe ~ ^ W *«-m =0 j = 1,NP ; 


4j+l 


P } - v*) a 4j^4 =0 j = l.NP 


< u3 to • vj) °4j+4 = 0 j = l ’ NP ’ 


(4.63) 


(4.64) 


and finally, for n = 4 + 4NP + i: 


P (u to ' V) ff 4+4NP+i 


= 0 i = 1, NG 


or, 


(u to ' v) a 4+4NP+i 


= 0 i = 1,NG 


(4.65) 


Now, substituting the values of and from Eqs. (4. 19), ( 4 . 20), (4.24), 
(4.25) and (4.28) into the generalized compatibility equation (4.15) yields: 


L = (p (Xj + pu a^) du + (pu Oj) dv + (a ^ + ua 4 ) dP + (uoj - a ucr 4 ) dp 
NP 

+ Z |<P J CT 4J+1 + ^ d “ 3 + & “ 3 0 4j + 3> dv> 4 ,u j « 4J+l ) dp) 


N ...... 

+ (p j u j a 4 . +4 ) dh j ] + 5^ pu <y 4+4NP+i d* L + 1^ a l +a 2 £ P 3 A 3 < u - u3 ) 

i = l 


NP 


j=l 


NP . NP NP r .J j 

+ a 3 !C p3 A3 < v “ v 3 ) + a 4 W»| ‘ £ P 3 a3 b {) + 2 I y V a 4j 
j = l j=l j=l 


(4.66) 


4j+l 
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- P j aJ < u - uJ > a 4j+2 ' p3 AJ (v ' vJ) a 4j+3 + P J 4 ff 4j+4] 


NG 

11 w i ®4+4NP+i 
i = l 


dx 


0 . 


(4.66) 

(Continued) 


Equation (4.66) is the expanded form of the generalized compatibility 
equation and is used to determine the compatibility equations that are valid 
along each characteristic. 


To determine the compatibility equations that are applied along gas 
streamlines recall that along gas streamlines 


*1 = I 
dx u 

Therefore, Eq. (4.57) becomes: 

£ + 0 

or 


Equation (4.59) becomes 

a 2 S' • a 3 + 4 = 0 

or 

dx u 

(Jq ~ 1 ri (T« — a • 

2 dy 3 v 3 

Along a gas streamline: 

£ 0 j = 1, NP . 


(4.41) 


(4.57) 


(4.67) 


(4.59) 

(4.68) 
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Therefore, Eq. (4.63) becomes: 

(u J ^ - v j ) a 4j+1 * 0 j = 1,NP (4.63) 

or 

a 4j*l * 0 J = 1,NP * (4.69) 

Similarly, Eq. (4.64) becomes 

(u j - v*) or 4j+4 =0 j = l.NP (4.64) 

or 

V 4j+4 = 0 j = 1 > NP * (4. 70) 

Recalling Eq.(4.69), Eq.(4.62) yields 

0 • A 

Si fa * <“ J dJ • yJ » ff 4j+3 * 0 J*'- NP (*•«> 

or 

a 4j+3 = 0 J * • (4.71) 

Similarly, Eq. (4.61) yields 

-0 

S^j+1 + < uJ to ' vJ) a 4j+2 = 0 j = 1,NP (*•«> 

or 

a 4j+2 = 0 J = l ' NP • (4.72) 

Substituting Eqs. (4.67)through (4.72) into Eq. (4.66) 

2 

L = o 3 du + pu a 3 dv + (^ a 3 + u a 4 ) dP + (-a u a 4 ) dp 
NG NP 

+ 1C P u 7 ‘ ^3 dx + 

i=l j=l 
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NP . . 

+ ^ p * A * (v - v^) ffj dx + 

j=l 


N£ 

a 4 ( *1 ' £ 


Bj) dx - ^ O', 

j=l i=l 


NC 


and rearranging terms fields 


| 2 N p f i 

L = do + pu dv + — ■ dP + p^ I ~ (u - u^) + (v - v^) I dx J 

I , NP . . . I 

+ |u d p - a u dp + (4<j - p* A^ B j) dxj cr^ 

l j=l I 

NC , 

+ Z)| a 4+4NP+i <P U - w. dx)j = 0 . 


4+4NP+1 dX = ° ’’ 


(4.73) 


Since the remaining multipliers, a, are arbitrary, their coefficients must 
be zero. Therefore, the coefficient of a ^ yields 


NP 


du + pu dv + ^ dP + £ p j A- 5 (u - u j ) + (v - v j ) j dx = 0 


or 


Since 


j=l 


pu du + pv dv + dP + p* A^ |(u - u^) + ~ ( v “ v^)| dx = 0 


j=l 


2 2 , 2 
q = u + v ; 


(4.74) 


or 


2q dq = 2u du + 2v dv 


q dq = u du + v dv 


(4.75) 


Substituting Eq.(4.75) into Eq.(4.74), and dividing through by p yields: 

NP 

q dq + ^p + p iC ^ [ (u ' + u (v ’ V ^] <** = 0 * 

j=l 


(4.76) 
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NP 


u dP - u dp + dx p^* b| dx = 

j=l 


= 0 


or 


»|> NP . , 

dP - a^ dp + •— dx - — 23 P ^ b| dx 

“ j=l 

and finally for the coefficient of ®4^4j$P+i 

f(pu dX - dx) = 0 


= 0 


or 


i=l 


pu dX - dx = 0 i = 1, NG 


(4.77) 


(4.78) 


Equations (4.76), (4. 77) and (4.78) are the compatibility equations that 
apply along the gas streamlines. 

To determine the compatibility equations that apply along Mach lines, 
recall that along each Mach line 


and, 


•j^ = tan (9 + a) 

djr y v 
dx * u * 


(4.56) 


Equation (4.60) may be rewritten in the form: 


-* 2 v * 0 


(4.60) 


Therefore, 


or 


CTj - a a 4 =0 


a l = a a 4 * 


(4.79) 


Equation (4.57) becomes 

°l to + °2 ,u to ' v) ’ 0 
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Using Eq. (4.32) 


-a dy/dx 
°2 “ (u dy/dx - v) 


S = u - v ; 
dx 


the expression for or ^ becomes 


= 


-a <y ^ dy/dx 


Equation (4. 58) becomes 


or 


o~ = 


• u *- v > < ’3 


a l a °4 


3 S 


Recall that along a Mach line 


J & . vi 


u- - v J / ° j = l.NP 


Therefore, Eq. (4.63) yields 


- A o.... =0 Jal.NP 


3 ' 


or 


’«i+i » j * i,np • 


Applying Eq. (4.82) to Eq. (4.6)) 

0 

6 ,1 +(uJ |J ' vJ)<y 4i+2 = 0 i =l * NP 


or, 



ff 4j+2 = 0 j = 1,NP • 


Simile ?ly, Eq. (4.62) yields 

r0 


* < ui S’ - vJ » °4j+3 = 0 j s ‘• NP 
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or. 


<r 4J+3 =0 j = 1,NP . 


Equation (4.64) yields 




or 


a 4j+4 = 0 


= 0 j = 1, NP 
j = 1, NP . 


Finally, Eq. (4. 65) yields 


or 


<“ & - »> '4<4NP« * 0 t = I - NG 


a 444NP+i = 0 i = 1,NG . 


(4.84) 


(4.64) 

(4.85) 


(4.65) 

(4.86) 


Substituting Eqs. (4.79) through (4.86) into Eq. (4.66) 

2 


„ 2 

2 pu a a 
L * (pa <J 4 = 


4 dv. , . P“ » , . / 

— s — dv + (- 


a" dy/dx 


uo 4 ) 


dP 


+ (« dp + ^ a 2 a 


a 2 a 4 dy/dx NP 


/ 


2. NP 


-g -£p j A j (u - J) 

j=l 


a a 


. / NP A) 

( v - v ^> + °4 Pi ‘ 2 P * B l)j 

j ^ 


dx = 0 ; 


and rearranging terms yields 

*• ■ |K - V £)*• ♦“£*♦(.- 4 *.)dp 

a 2 + T S'> i A> [ ,v ‘ vi) ' ( “' ui) + +1 P j A J B J ijdxjo 4 s 0 
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Since the multiplier, , is arbitrary, the coefficient must be zero. Therefore, 


(p * 2 - “r- &) d “ +e ¥- dv + - t &) 


dP + 4ev»_d* 

y 


+ j(v - \r*) - (u - u^) ^ Jdx + ipj dx - b| dx =0 (4.87) 

j=l J j=l 


The co^':'icier.t of du may be rewritten in the form 


_ 2 Pu a dy __ 2 
P* ' S' dx = Pa 


(S ^| tZj g ) = pa2 


(4.88) 


The coefficient of dP may be rewritten in the form 


2 2 2 2 
a dy uS - a dy/dx _ u~ dy/dx - uv - a dy/dx 

S dx S ~ S 


a 2 dy (u 2 - a 2 ) dy/dx - uv 
Sdx = S 


(4.89) 


Recalling that 


dy -uv +a 2 ^M 2 - 1 

dx 2 2 

a - u 


and rearranging terms 


. 2 2. dy , 2^/».2 , 

(a - u ) r* = - uv + a w M - 1 


(u 2 -a 2 )g- = uv+a^M 2 -! 


finally, 


(u 2 - a 2 ) ^ - uv = + a 2 ^M 2 - 1 . 


(4.90) 


c -X 


4-25 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC aR D867400-I 


Substituting Eq. (4.90) into Eq. (4.89) yields 


M 2 - 1 


= ±a. 2 Jf 

S dx S 


Substituting Eqs.(4.88) and (4.91) into Eq. (4.87) 


pv du + pu dv + ^M 2 - 1 dP 


S 


2 NP . 


+ ^_ dx 

y 

np . 


+ P* A^(v - v^) dx - (u - u^) dyj + «J»j dx - ^ f? B j dx 

j=l j=l 

2 

d iding through by a p/S and rearranging terms yields 
udv - vdu +WM“ - 1 — +(— + — Sdx 




NP 


+ p^ T(v - yj) dx - (u - u^) dy - b| dxj 


= o 


Recall that 


u = q cos0 , 
v = q sin 9 , 


and dy/dx = tan(9 + a) along each Mach line. 


Therefore, S dx may be rewritten in the form 

S dx = ( u - v) dx = (q cos9 tan (9 + a) -q sin9) dx 

and upon rearranging terms 

Sd* = qdx[co.9 “”(»+“) . sln9 £2£l*Ia.l 
cos(0+a) cos(0+a)J 


4- 26 


(4.91) 


; (4.92) 


(4.93) 

(4.46) 

(4.47) 

(4.56) 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC TR D867400-I 


S dx = 


_ad£_ I 
(0+a) ' 


co? 


cos0 sin (9 + a) - sin 9 cos 


(e+d)J 


(4.94) 


To simplify the above equation consider the trigonometric identity: 

sin|(9 +a) _+ sj = cos0 sin(6 + a) + sin0 cos(0 + a) . (4.95) 


Therefore, 


S dx s 


— — sin 1(9 +al - ol = q dx t q) > 

cos(0+a) * * cos (9 +a) 


or 


„ , — o sma , 

S dx = + — dx 


(4.96) 


co8(0 +a) 

From Eqs. (4.46) and (4.47) above the following expressions can be written: 


and 


du s - q sin0 d0 + cos0 dq 


dv = q cos0 d0 + sin0 dq 


(4.97) 

(4.98) 


Therefore, the first term in Eq. (4.93) may be written as 
2 2 2 2 

u dv - v du = q cos 0 d0 + q dq sinQ cos0 4 q sin 0 d0 - q dq sin0 cos0 


or 


and finally, 


2 2 2 

u dv - v du = q d6 (cos 9 + sin 0) 


u dv - v du = q d6 . 


(4.99) 


Substituting Eqs. (4.96), (4.99) and (4.50) into Eq. (4.93) yields the final form 
of the compatibility equations that apply along Mach lines 


q^ d0 + cota ~ + 
NP 


( 60 3 in 9 + ^1 \ h 3 sing \ 
^ pa / \ co8(0 +ay 


dx 


(4.100) 


+ ^ V p 4 |(v - v^) dx - (u - u^)|dy - B-j dx j = 0 

j=l * 
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Finally, to determine the compatibility equations that apply along the particle 
streamlines recall that along each particle streamline 

^ to ' ^ = 0 » 


or 


and that. 


. v" - 


dx i 
u J 


= tan© 


j 


(4.42) 


±L A Z 
dx ' u 


We may immediately rewrite Eqs. (4.79), (4.80) and (4.81) which were derived 
for Mach lines. 


Equation (4.79) 



Oj = a a 4 , 

(4.101) 

Equation (4.80) 

a 2 a. , 

a 2 S dx * 

(4.102) 

Equation (4.81) 

a 2 o 4 



a 3 ~ s 

(4.103) 


Then applying Eq. (4.42) to Eq. (4.61) yields 


iK j+1 WJ - 



0 
4j+2 


=0 j = l.NP , 


(4.61) 


or 


UK*! J * >' NP 


°4j+l * 0 j = 


(4.104) 


Similarly, Eq. (4.65) yields 

(-M 


a 444NP+i = 0 1 = l.NG 


(4.65) 
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or, 


a 4+4NP+i = 0 1 = 1,NG * 


Substituting Eqs. (4.101) through (4.105) into Eq. (4.66) 

L = fpa‘o. -fia-aia 4 ^du+fi^-<T J ,dv+/ua > , -^r^-^irJdP 


= a* 




'4 S w 4 dx/ u “ T ~ s - u 4 “ v ' \ 4 S 

NP 


+ (u a a 4 - u a o 4 ) dp I P J u J j +2 du J + P J u J <^ 4 j +3 dv 1 

2 2 NP 

. ; .i 6pv a a o A , . 

+ pJ u J dh J j + dx g — dxS^ aJ (u * uJ ) 

j=l 

2 NP NP 

a (J 4 * ■ ^ ill 

+ — g 2_J p* A J (v - yj) dx + a 4 ‘I'j dx - a 4 /P J A J B j dx 

j=l j=l 

NP 

+ S[' pJAJ( “‘“ J)a 4j+2' piAi(V " vj,0 4j+3 +pi ' , i Cr 4j+4] dx 1 C 

and rearranging terms yields: 

L ■ J <» 2 (* -Sf) da+ p a 2 l d - + (“-T^) dP + 6 - £ r' <ta 

2 NP 

+ ^ 2 P^ A ^ j(v - v^) dx - (u - u*) dy - B j dxj + ijjj dx 

j=l a 


NP 


+ 2j[p j u^ du^ - p^ A^ (u - u*> dxl a 4 j + 2 

j=l J 

NP 

+ 2^ £p^ U^ dv^ - p^ A^ (v - v*) dxj <^ 4 j +3 


j=l 

NP 


+ [V dh^ + pj ^ dxj o A m - 0 


j=l 
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Since the remaining multipliers, cr, are arbitrary their coefficients must be 
zero. Therefore, the coefficient of yields 

Pa 2 (> -|S)' i " + Pa 2 | dv 4 * 4 &)«* +i£ r !d * 


NP 


+ p^ A^J(v - v^) dx - (u - u^) dy - b| dxl + <|»j dx 


= 0; (4. 


the coefficient of ff 4j^2 

NP 


or 


W * » t • • • • * 1 

y* lp^ u^ du^ - p** (u - u^) dx I = 0 ; 

j=l 1 

u^ du^ = A^ (u - u^) dx j = 1 , N P . 


(4, 


Similarly, for the coefficient of 


NP _ ... 

^ J"p^ u^ d\P - p^ A^ (v - \r*) dxl = 0 
3=1 J 


or 


u^dvj = A j (v-\r*)dx j = 1,NP ; 


(4. 


and finally for the coefficient of ^ 

NP 


or 


E[p j u^ dh^ + p^ ^ dxj = 0 
u^ dh^ = - 4*9 dx j = 1, NP . 


(4 


Recalling that 


.4 

= | A* (T j - T) + -py (T j ) - J T 4 J ; (4, 
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Eq. (4.110) becomes 
u^ dh^ = - 


■| C j (T* - T) + y j ji? J (T^) - a* T 4 || 


dx j s 1, NP (4.111) 


Before evaluating Eq. (4.107) it is necessary to examine Eqs.(4.57) through 
(4.60) with S = u dy/dx - v non -vanishing (S ^ 0) as in the case along a particle 
streamline. 


Equation (4.60): 


or 


S Oj - a S o 4 = 0 


Oj = a o 4 ; 


(4.112) 


Equation (4.58): 


or 


o, = S <7, 


= tf/S 


finally. 


2 

a o t 


<7, 


(4.113) 


Equation (4.57): 


or 


= 0 
2 

<7. i a o A j 

1 dy 4_ dy 

°Z = " S dx = ' S dx 


(4-114) 


and finally 

Equation (4.59): 


<7- 4 ^ - O-, + S a a - *J . 
2 dx 3 4 


(4.59) 
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Substituting Eqs. (4.112) through (4.114) into Eq. (4.59) fields 

® 2 °4 JdY\ ^ ”4 
‘ ~ S - (to) ‘ “ S - + S ®4 0 


or 


[- 


* 2 ($-» 2 ^ 2 


I-*' 


In the above equation, either is zero or its coefficient is zero. 
Along particle streamlines: 

dy _ yJ 

dx - j 
u J 


Assuming cr^ = 0, Eq. (4.115) becomes 

.2 


1 - 1 ?) 


Recalling that 


or 


S = u£-v 


S = u — r - v / 0 ; 
uJ 


Eq. (4.115) may again be rewritten in the form: 

2 

2(v^)f 2 , 2/v^\ , v^ . 2 n 

- a - a + u I — r I - 2 uv —r + v =0 


VuV uJ 


or 


$ 


<u 2 -a 2 >(4) - 2 uv 4 + (v 2 - a 2 ) = 0 

u^ 


Solving for — r- , 
u^ 


J 2 uv + ^4 u 2 v 2 - 4 (u 2 - a 2 ) (v 2 - a 2 ) 


u J 


2 (u 2 - a 2 ) 
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upon rearranging terms, 


* 


. 2 2 2 22 2,2 2 4 

uv + Vu v - u v +a v +a u -a 


u' 


2 2 
u - a 



This equation has the same solution as Eq. (4.43), that is: 

v> 

— r = tan (9 + a) . 
u^ 

However, for a particle streamline, 

v j j 

—7- = tan9 J ; 

u^ 

therefore the assumption that the coefficient of a ^ in Eq. (4.115) vanishes 
does not apply along particle streamlines, hence the multiplier must be 
equal to zero and Eq. (4.107) is not valid along particle streamlines. In 
summary, the compatibility equations that apply along particle streamlines 
are: 

. it = tane^ j = l.NP 
* U J 

u^ du^ = A* (u - u^) dx j = 1,NP (4.108) 

u^ dv^ = (v - v^) dx j = 1,NP (4.109) 
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and 


dh^ = - 


j A J C J (T J - T) + 


-J-T [^(T j ) 4 - a j T 4 ] 
m J r J L J < 


dx j = 1, NP (4.111) 


4.2.2 Enthalpy -Entropy -Velocity Form of the Compatibility Equations 

Recall the pressure-density-velocity form of the compatibility equations 
that apply along gas streamlines 


NP 


q dq + -p- + P^ J(u - u J ) + J (v - v J )j dx = 0 
j=l 

tjj NP 

dP - a 2 dp + dx - £ ^ P^ A * B j dx = 0 

“ j=l 


(4.76) 


(4.77) 


pu dX. - w^ dx = 0 i = 1,NG (4.78) 

Equation (4.76) may be written in the enthalpy -entropy form using the thermo- 
dynamic relations 

NG 

T dS = dh - dX. (3.88) 

1 = l 

2 

H = h + ^- (3.68) 

and the definition of the local speed of sound. 


From Ref. 9 it is shown that since the particle velocity and temperature 
do not change through a discontinuity and since there is a finite relaxation 
time associated with changes in the particle properties, then an infinitesimal 
weak disturbance must travel at the gas-sonic speed. That is to say, the local 
speed of sound of the gas-particle mixture is unaffected by the particle and is 
identical with the speed of sound of the gas; hence, it can be expressed by the 
following 

a 2 --- yRT . 
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From Eq. (3.68) we may write 

dH = dh + q dq 
or 

dh = dH - q dq 


( * 6 ) 


Rearranging Eq. (3.88) 


dP 

P 


NG 

dh - Tds - jtL dX. , 
i - 1 


substituting for dh from Eq. (4.116) 


dP 

P 


NG 


dH - q dq - TdS - dX. 


i = l 


and rearranging terms yields 

dp 

~ + q dq = dH - TdS - ^ . (4.117) 

irl 

Substituting this result back into Eq. (4.76) yields the final form of the com- 
patibility equation 

NG NP 

dH - TdS - PjdX. + ^ £ pj [(u - u J ) + - (v - vM dx =0 (4.118) 

i=l j=l 

Equation (4.118) may be rewritten using the definition of th: local speed of 
sound. 


Solving for temperature 


2 . 2 
q sin Q 

yR 


(4.119) 
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and substituting the result back into Eq. (4.118) yields 
2 NG NP 

dH - q 25 ^- dS -£ M i d^i + £ X)P j A j [(u-u j ) +^(v-v j )]dx = 

i=l j=l 

Next, we examine Eq. (4.78). 

Equation (4.78) may be rewritten in the form 


or 


p q cos0 dX = w^ dx 
w^ dx 

^i “ pqcosO * 


i = 1, NG 
1, NG 


By letting 


Eq. (4.122) becomes 


or finally, 


w i 

X = — - 
i P ’ 


X . dx 

d *i 5 q'cosij i = 1 ’ NG 


qco.edXj = X i dx i = l,NG 


Equation (4.77) may be written in the enthalpy -entropy form using the 
thermodynamic relations 


and 


Rearranging Eq. (3.97) 


T dS 


+1 


dP 



|/dP _aidfi\ 

V R I 

‘ P P / 

1 

NG 

V ,, • 

cjx~- 

i 2- ^i w i • 

2 dp = 

P T dS 

(Cp/R - 1) * 
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and substituting fcr T from Eq. (4.119) yieldt 

dP - a 2 do - 

dF dp ' r (C - R) * 

r 

Replacing the first two terms of Eq. (4.77) with Eq. ( A .125) 

NP 

p S 7 r 8i "\ -T S Adx-lVpiA 1 B J . dx = 0 
y (C - R) a u l—J r 1 

p j=i 

substituting for from Eq. (3.104) 


2 2 NG NP 

. 1 V'pj A j B j ^ . o 

y(C - R) (C -R)u i l 

P P i=i j=l 


^nd using the relation X. = w\/p results in 


NG 


NP 


0 7 1NU IN Jr 

P " 3 ~, f ^~ n - g r s - * T r - ~ £ \ ax.iVpi A iBi dx 

y (C - R) (C - R) q cos0 / * i * u L-d^ 1 

P P i=l j=l 

Dividing through by pR/\Cp - R), and recalling Eq. (4.119) yields 
NG (C - R)NP 

Tds + i-^Fe2^ M i ; 'i dx -^lR— 2 piAi B l dx = 0 

; - 1 i_i 


j=i 


(4.125) 


(4.126) 


0 . 


(4.127) 


where 

B = c ~~^ R k * Aq j - ^ • A? j + f C j (T j - T) + . 3a r ■ W (T j ) 4 - Ct j T 4 ]! (4 . 

p ( A J m J r J *- J * 


103) 


The final form of the compatibility equations that apply along gas streamlines 
are summarized as follows: 


Equation (4.127) 
1 


NG 


T dS + 


q cosO 


£ p i*i 


dx - 


i = l 


< C P - R ) 

upR 


NP 


pj B? dx = 0, 


(4.128) 


j=l 
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Equation (4.118) 


NG 

NP 



dH - T dS y] 

+ iT* t J A j f(u - J) + J (v - 

)| dx = 0 

J 

(4.129) 

i=l 

M 1 




and 

Equation (4.123) 

q cos0 dX { = X. dx i = 1,NG (4.130) 


The pressure-density -velocity form of the compatibility equations that 
apply along Mach lines are 

q 2 as + coto ® + (ia-isaS + -)dx 

P \ y pa 2 '' cos(0 + a)' 

N p 

+ a' [(. - v') dx - (u - u') dy - B j dxj s 0. (4.100) 

j=l * 


Substituting Eq. (4.117) 


NG 


= dH - TdS -qdq - £ ^d*. 

i=l 


(4.117) 


into Eq. (4.100) 
q 2 d9 “ 


NG 

+ cota (dH-T dS - q dq - V fJL. dX) I 3 . ** S £L *3_2H)§ 

Ul 1 1 COS (»+«)' y 


NP 


+ ^) dx+ pE piAj 

pa j=i 


(v - v*) dx - (u - u^) dy + ^ 


8 in a 


cos(8+a)q sin 


— = B\ dxl = 0 ; 
in a J 


4-38 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC TR D867400-I 


dividing through by q 


2 


and expanding wherever possible yields: 


t cota JIT . cota . cota . cota 

d9 + — r- dH + — T aS + dq + — =- 

q q H q 


NG 




i = l 


— 6 sina sin9 , — 
+ dx + 


sina 


cos(9 + a) y qpa cos(9+a) 


1 n . 1 

- dx + — 


pq 


,P J A J 


j=l 


(v - v^) dx - (u - u 3 ) dy 


Bj dx 


q sina cos (9 +a) 


= 0 


(4.131) 


Equation (4.131) may be rewritten in the following manner: 


Recall that 


and 


2 . 2 

T = a ,-si q q 




y R 
NG 

cT /5 ~ - i 5>i*i 

i = l 



Therefore, the third term of Eq. (4.131) may be rewritten as 


or 


£2t|L TdS = cotaq 2 sin 2 a dg 

q q yR 


cota sina co sa , c 

1 X do - do 


y R 


The seventh term may be rewritten as 


sina dx 

2 — 
q pa cos(0 + a) 


r i . l 

SinCt (cT7RT7j Z-f M i W i 
l p i=l J 


dx 


3 2 — 

pq sin a cos(9 +a) 


(4.119) 

(4.104) 


(4.132) 


(4.133) 
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sina 4*2 dx 


dx 

C_/R 


2 — ~ 3 — 

qpa cos(0+a) q sina cos (0 4 a) 


(A. 133) 


A portion of the eighth term may be rewritten as 


(v - v^) dx - (u - u*) dy = 


= dx |(v - - (u - u^) ^ J 

= dx j(v-v^) - (u-u^) tan(0 4a)| 

= dx f(v - yj) c ? 8 Ji.+ a ). - (u . u j) 8 i S(g +°i | 
| cos( 0 +a) cos( 0 +a)j 


finally, 

(v - v^) dx - (u - u^) dy = — f(v-v^) cos(0 4a) - (u-u^) sin(0+a)j . (4.134) 

cos(0 +a) * * 

Substituting Eqs. (4.132), (4.133) ana (4.134) back into Eq. (4.131 yields 
the final form of the compatibility equations that apply along Mach lines 

. cota , . sina cosa dS — cota dH — 6 sin9 sina dx 

<10 ±— dq 4 -p + 2 + = 

q ya ycos(9+a) 


- _ 2 


pq cos 


NP r i 

— — — p* A? + (v - v^) cos (0 +a) + (u - tr) si n(0 +a) + — 

(0+a) . | 

J“ x Mr. 


« ul 


iivjr 


l_ i=i 

3 — 

q sina cos(9 4a) 


(4.135) 
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4.2.3 Compatibility Equations for Equilibrium and/or Frozen Chemistry- 

Recall the enthalpy-entropy -velocity form of the compatibility equations. 


Along gas streamlines: 

NG (C - R) NP 

T as ZjMi - -fp*- £> A * »{* ■• 

and 

NG NP 

dH - T dS ^ +^£ pj A3 [ (u_uj) + u (v " vj) ] ^ = °* 

i=l j=l 


(4.128) 


(4.129) 


Along Mach lines: 


, cota , . sina cosa dS — cota dH — 6 sin9 sina dx 

d9 ±— — dq + -5 + j + r 

q 7 q y cos(0 +a) 


, NP . r _ . _ B j . 1 

+ — 5 Z — 2j p* A J | + (v - \P) cos(e +a) + (u- u J ) sin(8 +a) + — -r — I 

Pq cosO+a)" [ q j 


NG 


, cota r - ' r 

± 2 iJ dX i + g— 


NG 

c^tn E Mi - 0 

_ P i = l 


(4.135) 


i=l 


q sina cos(9 +a) 


For any closed system which is in the state of complete equilibrium the 
terms NG 

YV dx. 

k 1 1 

and NG 

Em> 

i = 1 

must be zero (Ref. 40). Furthermore, for chemically frozen flowX. and dX. 
are obviously zero since there are no changes in the chemical species con- 
centration. Therefore, for the case of equilibrium and/or frozen chemistry 
the above equations reduce to the form: 
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Along gas streamline 8: 


(C - R) 

TdS --fe- 


NP 


and 


]^p j A j B J ‘ dx = 0 , 
j=l 


NP 


dH - T dS + [(u - J) + ~ (v - v*)j dx 

j=l 


= 0 . 


(4.128a) 


(4.129a) 


Along Mach lines: 


.. . cota , . sina cosa dS — cotadH — 6 sm0 sma dx 

d9 + — — aq + —5 + 2 *•“ _ 

^ o~ ycos(0 + a) 


NP * _ _ B j 

+ — =-^ — — p^ aJ |± (v-v^) cos(0+o) + (u-u^) sin(0+a) + \ 

pq 2 cos(0+o)" I 4 


(4.135a) 


4.2.4 Summary of the Compatibility Equations for Gas-Particle Flows 

The characteristics of the flow have been found to be the gas streamlines, 
gas Mach lines and particle streamlines (one for each particle species). 


Pressure-Density-Velocity Form (for Chemical Non -Equilibrium and Transi- 
tion Flow Applications) 

The slope of the gas streamline, 0, is given by 

^ = tan0 (4.41) 

and the compatibility equations which apply along gas streamlines are: 

NP 

q dq + ^ + 2 P J A j [(« -u j ) + £ (v - v j )l dx = 0 , (4.76) 

j=l 1 J 

<l> NP 

dP - a 2 dp + dx - £ £p j A j B] dx = 0 (4.77) 

j=l 
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and 

pu dX - dx = 0 i = 1,NG (4.78) 

The elope of the Mach lines (left running characteristics and right 
i inning characteristics) is given by 

^ = tan (6 +a) -56) 


and the compatibility equations which apply along each Mach line are: 


d0 + cota 


dx 


dP t lb q sinQ ^ ^1 \ /y q sina \ 

P \ ^ pa /\ cos(0+a)/ 

NP . 

|(v - v^) dx - (u - u*) dy - B j dxj = 0 


(4.100) 


Equation (4.100) may be written in the more convenient form 

NP 

dP - 


d0 + cota 


Pq 


- 6»m8 sinadx ± . dx _ A> [+ (v - v*) cos(9 

y cos (0 +a) pq^ cos (9 +a) 4^ | 


+a) 


+ (u - u J ) sin(0 +a) + 


B 1 


. NG 

mEMi 

+ Lzi = 0 


q sina | q 3 sina cos(0 +a) 


(4.136) 


The particle streamline direction, 8^, is given by 


4* = ^ = tan0 j j = 1, NP (4.42) 

u J 

and the compatibility equations which apply along particle streamlines are: 

u j du 3 = A j (u-u J )dx j = 1,NP (4.108) 

u 3 dv^ = (v-v^) dx j = 1,NP (4.109) 
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and 


dh^ = 


fJc 


i (T^ - T) + — yy (T j ) - Ot* T 4 jJ 


dx j = l.NP (4.111) 


Enthalpy-Entropy-Velocity Form (For Chemical Equilibrium and/or Frozen 
Flow Applications) 

The slope of the gas streamline, 0, is given by 


& . 


dx 


= tanO 


(4.41) 


and the compatibility equations which apply along gas streamlines are: 
dH - T dS + p* A* [(u - u*) + ~ (v - v j )] dx = 0 

j=l 

(C - R) NP . . . 

T dS - -^5ir 2> J aJ B 1 = 0 

j=i 


(4.129a) 


(4.128a) 


where 


B 


1 ■ cybi A? + 1 ci (T i - T) * ( T i, 4 - ai t 4 ]| 


(3.103) 


and 



Uj L . - 

m' 1 (r^) 


c j = tsi 

v i> 


The slope of the Mach lines is given by 

^ = tan(0+a) 


(3.46) 


(3.83) 


(4.56) 
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and the compatibility equations which apply along each Mach line are: 


. cota . sina co»a dS — cota dH — 6 s:n0 sina dx 

do + — “ — dq + 5 + s + — 

** ^ q y co8(9 +a) 


dx 

pq 2 cos 


NP f 

- £ P^ | + 

(9 + a) ^ [ 


(v - v^) cos(0 +a) + (u - u^) 8in(0 +a) + 


4-1 

q sina I 


(4.135a) 


The independent variables are x, y, q, 0, P, p, u^, v^, and p^ in the 
Pres sure -Density -Velocity form, and x, y , q, 9, H, S, u^, v^, h^ and p^ in the 
Enthalpy-Entropy-Velocity form. 6+4NP unknowns are required to com- 
pletely define the gas -particle flow at a given location in the flow field. The 
above equations provide only 6+3NP compatibility relations. This is because 
a compatibility relation does not exist for the particle density and results 
from the assumption that the particles do not contribute to the pressure acting 
on a control volume in the gas -particle system. As originally shown by Kliegel, 
(Ref. 10) this can be resolved by utilizing the following incompressible flow re- 
lation to obtain the particle density. 


. • * { * • 6 * * • 6 *| 
dm^ = (2jr) 6 (y^) dy^ - (y^) dx^j 


(4.137) 


and 

6 takes on the values 

6 = 0 for two-dimensional flow 

1 for axisymmetric flow 

4.3 FINITE DIFFERENCE SOLUTION OF THE COMPATIBILITY EQUATIONS 

To solve the system of compatibility equations it is first necessary to 
write these equations in finite difference form. Notice that the equations de- 
scribing the particle properties were derived for a discrete particle. To 
account for more than one particle size or species, n discrete particles are 
allowed in the numerical solution. The streamline of each additional particle 
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is an additional characteristic curve; and the solution is obtained by applying 
the particle compatibility equations independently for each particle along its 
streamline. 


The equations are cast in a form suitable for the calculation of an 
interior point in the flow Held. For special points, such as along the nozzle 
axis, the nozzle wall, the limiting particle streamline and a free boundary, 
certain boundary conditions in the flow field are known which allow some of 
the equations to be discarded. The characteristic net for an interior point 
is illustrated by Fig. 4- 2. 

The difference equations for the gas -particle system of compatibility 
equations in coefficient form is given below. 

The slope of the gas streamline, 0, is given by 


|*5-3 = tan'b 5 3 

and the compatibility equations which apply along gas streamlines are 


Pressure-Density-Velocity Form for Finite Rate Chemistry 

AP, 


l 5, 3 Aq 5-3 + 1T 


5-3 


5,3 


+ C 3p = 0 


and, 


AP 5-3 " a 5, 3 ^5-3 + C 2 p = 0 


*>5, 3 u 5, 3 AX i 5 _ 3 - \ 3 ^5-3 = 0 1 : 


(4.138) 


(4.139) 

(4.140) 

(4.141) 


where the barred values are averages over the step length, and the coefficients 
are defined as follows » 
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NP 

= - 1 £ 3 A 5,3 [ (u 5,3* u 5, 3* COS ®5, 3 + (v 5, 3 ' v 5, 3 )sUl ®5, 3] AL 5-3 


'P P 


5,3 j=l 


and 


1 


c - 1 

P 5, 3 


L "5,3 Ji-1 


NG NP 

£ \s\ 3 '2^5. 3 *5, 3 B i 


j=l 


5,3 


AL 


5-3 


^ 5,3 


AL is the step size along the streamline. 

Enthalpy -Entropy -Velocity Form for Equilibrium and/or Frozen Chemistry 


AH 5-3 = T 5, 3 AS 5-3 " C 3 


(4.142) 


and 


AS 5-3 = C 2 


(4.143) 


where 


C 2 ' 


C - U rNP 

-¥ )z 

5 - 3 / L J * 1 


aJ gJ 
P 5,3 A 5, 3 


5,3 


AL 


5-3 


P 5, 3 T 5, 3 q 5, 3 


The slope of the Mach lines (left running denoted by subscript 2 and 
right running denoted by subscript 1) is given by 


1, 2 = tanp^ 


(4.144) 


where 


P l ,2 = 0 1,2 + X 1 , 2 ° 1,2 
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and X takes on the values 


X 


1 


X 


2 


1 for interior solution 
1 for lower boundary solution 

G for upper boundary solution 

1 for interior solution 

0 for lower boundary solution 

1 for upper boundary solution. 


The compatibility equations which apply along each Mach line are: 


Pressure-Density- Velocity Form 


A 0 1,2 + Q 1,2 AP 1,2 G l , 2 - C 1 l ,2 = ° 


(4.145) 


where the coefficients are defined as follows 

cosot 


Q 


1,2 


I4.2 


' 1,2 


Sm0l l, 2 q l, 2 P l, 2 


6sin6 1 2 sina 1 2 Ax t 2 

7i t2 c °8Pi >2 


and 

Cl 


1,2 


NP . r 

X) p 3 i. 2 A i, 2 [± 

j=i 1 


jki_i 

«l,2 ,ina l,2i 


<v l,2" , 1.2 ) co ‘Pl,2 +, "l,2-“i,2' sln Pl,2 


Ax 


~ — r 2 

p 


L± 


NG _ 

£ % 2 V i, 

i=l l ' c 1( 


Ax 


1^2 


1, 2 q l, 2 cosp l, 2 q l, 2 8ina l, 2 COSp 1, 2 C p, , 

1 


'1,2 
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Enthalpy -Entropy- Velocity Form 


A0 1,2- Q 1,2 Aq l,2- B l,2 +CI 1,2 AH 1,2 +G 1,2- C1 1,2 = 0 


where the coefficients are defined as follows 

cosa 


Q 


1,2 


LA 


sin0t l, 2 q l, 2 


B 


1,2 


sina 1 2 cosa 1 2 AS t 2 

a l,2 *1, 2 


cosa 


Cl 


LA 


12 — —2 

sino l,2 q l,2 


and 

Cl 


G , 6 sin9 1 2 sina 1 2 Ax t 2 


1,2 


y l, 2 cos Pi, 2 


1,2 


NP 

XX 2*1, 

j=i 


- < T 1, 2 ' v l, 2> co, f 1, 2 + '“1, 2 • u l, 2> 8 ‘" 5 1. 

D 1 2 

+ — LA 

q l,2 8in0t l,2 


Ax 


LA 


f‘l.2' l l.2 cosf l,Z 


Finally , the particle streamline direction, 0^, is given by 


X 4-3 = tanS? 
Ax 


4-3 = tanOi , j = 1, NP 

J 4,3 J 
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and the compatibility equations which apply along particle streamlines are; 


Au 4-3 

• «i. 3 

AxJ 

4-3 

j = 1, NP 

(4.148) 

< 

1 

* ^4,3 

^ 4-3 

i 1, NP 

(4.149) 

*>4-3 

* a i,3 

^i -3 

j = 1. NP 

(4.150) 


and 


where the coefficients are defined as iollows 

C2 j . A j lu 4,3' a 4, 3 1 
“• 3 4 ' 3 < 3 

j 


C3i . . A 1 


4.3 


4,3 


4,3 


and 


C4.j 

C4 4,3 


(4.150a) 


3 A 4, 3 C 4, 3 (T 4, 3 ~ T 4. 3* + “J [^4, 3 (T 4, 3> ” a 4 f 3 T 4, 3M 


4.3 4^3 


u 


4,3 


The difference form of Kliegel’s incompressible, flow relation is 

*“l, 2 = ,r •'* P \, 2 (“l, 2 <*1, 2 :> M. 2 ‘ 7 1,2 ^i, 2 > 6 ** 1, 2] • ,4 ' 151) 


An iterative solution is employed to determine the properties of the flow 
at the new point. For the first pass of the solution the Narred /alues are ap- 
proximated by the conditions at the known points, e.g., 


9 


1 " 


e 


1 


After the appropriate set of equations has been solved a new estimate of the 
barred values is made, e.g., 0 _ 
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The iterative process is continued until the desired convergence is 
reached . 

The mechanics of the numerical solution for a typical point are quite 
involved. Section 7 presents the mechanics of the numerical solution for the 

different types of points encountered in a typical gas -particle flow problem. 

4.4 DEVELOPMENT OF THE PARTICLE DENSITY RELATIONSHIPS 

Kliegel's incompressible flow relation 

• .6 . *5 .1 

u^ (y^) dy* - v^ (y^) dx^j (4.137) 

will now be used to develop the particle density relationships required to 
solve for the particle density at each of the different types cf points encountered 
in a typical gas -particle flow problem. The particle density relationships will 
be derived for an inferior point, a lower boundary point, a particle limiting 
streamline point and an upper boundary point. The particle density relation- 
ships for the remaining type points (such as a shock point) will not be developed 
because each is subject to the same boundary conditions as one of the above 
four cases and therefore yield the same results. 

4.4.1 Interior Point Solution 

Consider the characteristic net for an interior point 


dm^ = (2y) 5 


Normal 



y / V Left- Running 
Characteristic Line 



Right-Running 
Characteristic Line 
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Applying the principle of conservation of mass to the above system, we may 
write 


m l_2 ~ m l-3 + m 3-2 


(4.152) 


Each term in the above equation may be solved for by applying Eq. (4.137) 
between the proper points. 


Points 1-3: 


m, 


Points 3-2: 


m 3-2 = 


and point 8 1-2: 


m l-2 = 


f^i u i +p 3 u 3 ) (y^t - y 3 y 3> 

(PjVi + P 3 v 3 ) (y j + y 3 ) 

2 2 6 

2 2 6 

4 f(P 2 o 2 +Pj» 3 ) (y 3 y 3 - y 2 y 2 ) 

#>2*2 + P 3 V 3> 'A * A 

l 2 2 6 

2 2 6 

2: 

>(P 1 U 1 +P 2 U 2 ) (y l y l ‘ y 2 y 2 } 

(Pl v l + P 2 V 2 } (y6 l +y 2 ) 

2 2 6 

2 2 4 


( x 3 " x 2^ I ( 4 * 154 * 


The particle density relationship at point 3 may now be determined by 
substituting Eqs. (4.153), (4.154) and (4.155) into Eq. (4.152) and solving for p y 
Equation (4.152) becomes: 

^l u l + P 2 U 2> {y l y l ' y 2 y 2> ' (P 1 V 1 +P 2 V 2 ) (y l + y 2 } (x l * x 2> 

= +p 3 « 3 ) (y t yj - y 3 y 3 ) - (Pi v i + p 3 v 3 ) (yj + y 3 ) ( x i - x 3 > 

+ (P 2 u 2 +P 3 u 3 ) (y 3 y 3 - y 2 y 2 ) - (P 2 v 2 + p 3 v 3 ) (y* + y 3 > (x 3 - x 2 ) . (4.156) 
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To solve for p 3 : 

expand Eq. (4.156) 

p l u l (y l y l • y 2 y 2> ‘ P l v l (y l + y 2> (X 1 • x 2 ) +P 2 U 2 (y l y l ' y 2 y 2^ 

” P 2 V 2 (y I + y 2> {x l- x 2 ] = P 3 U 3 (y l y l ' *3*3* 

- P x vj (yj + yj) (xj - x 3 ) +p 3 u 3 (y 3 y 3 - y 2 y 2 ) +P 1 « 1 (yjyJ - y 3 y 3 ) 

+ P 2 U 2 (y 3 y 3 * y 2 y 2 ) * P 2 V 2 {y 2 + y 3 } (x 3 * x 2 ) * P 3 V 3 (y l + y 3> (x I ~ x 3 ] 

- P 3 v 3 (y£ + y 3 ) (* 3 - x 2 ) •» 

combine terms to form coefficients of P^.P^ and p 3 

P l[ U l (y l y l - y 2 y 2> * V 1 (y l + y 2> (X 1 ' X 2>] +P zl u 2 (y l y l - y 2 y 2> 

- v 2 (y l + y 2> (X 1 ' X 2 ) J = P 3 | u 3 (y l y l ‘ y 3 y 3 + y 3 y 3 * y 2 y 2> 

- v 3 (yj + y 3 ) (Xj - x 3 ) - v 3 (y* + y 3 ) (x 3 - x 2 )j 

+ P Y [u^y^t - y 3 y 3 ) - v^yf + y*> (Xj - x 3 )| + P 2 [u 2 (y 3 y* - y 2 y* } 

" v 2^ y 2 + y 3^ ^ x 3 ‘ x 2^| ’ 

isolate p ^ and its coefficient 

P 3 u 3 (y l y l • y 2 y 2> " v 3 [ (y 1 + y 3> (X 1 • X 3 ) + ,y 2 + y 3> {x 3 * x 2>] } 

= pj [u^y? - y 2 y 2> - v i (y i + y 2> (x i _ x 2> * u i (y i y i - y 3 y 3> 
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+ v i (y l +y 3> (k 1 - x 3 )] +P 2 [u 2 (y 1 y^ - Y^) - Wyf + y*> (Xj - x 2 ) 
- u 2 (y 3 /3 - y 2 y 2 ) + >^ 2 (y 2 + y 3 ) (x 3 - x 2 )] ; 

combine ; orms on the right hand side of the equal sign 

= ‘ y 2 y 2* ' v l|^ r l + y Z l ^ X 1 ' X 2 ) ‘ * y l + y 3* (X 1 ' X 3^|| 

+ P 2[ U 2 (y l y l * y 3 y 3 } * V 2 l (y l + y 2> (C 1 ' X 2 ) - (y 2 + y 3> (x 3 ' VjJ = 


and finally solving for 


P 3 = . 6 6, I. 6 . 6, . . . , 5 . 6 W ,1 I P 1 

u 3 y l y l " y 2 y 2* ' V 3l (y l +y 3 ) (X 1 ” x 3 + (y 2 + y 3 ) { 3 “ x 2 ) | 1 1 

- v,[(y* + y 2 ) (Xj -x 2 ) - (yj + y*) (Xj - x 3 >]| + P 2 j U 2 {y l y l * y 3 y 3> 


, 6 6 . 
1 ( y 3 v 3 - y 2 y 2 ) 


- v zl (y l + y 2> (x l • x 2 ) - (y 2 + y 3> (X 3 • X 2 ) ] |] ‘ 


(4.157) 


4.4.2 Lower Boundary Point Solution 

The lov er wall point and free boundary point are the two types of lower 
boundary points encountered in a typical gas -particle flow problem. To derive 
a relationship for the particle density at a lower boundary point it will be 
assumed that: (1) the lower boundary runs parallel to the x-axis (for two-phase 
only), and (2) the particle streamlines run parallel to the lower boundary when 
in the vicinity of the lower boundary. These assumptions are valid in nearly 
all of the flow problems encountered. 
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Consider the characteristic ne. for a point on a lower boundary 



Applying the principle of conservation of mass to the above system, we may 
write 

nij g = nij j tnij j . (4.158) 

Furthermore, from the above assumptions m^ v^ and v^ may be set equal 
to zero. 


Equation (4.158) then becomes 

m^ 5 = ”*1-3 * (4.159) 

Each term in the above equation may be solved for by applying Eq. (4.137) be- 
tween the proper points. 


Points 1-5: 


m l-5 = 


[<p u tp u ) , 1 

l T p l v l X 1 “*5 J 


(4.160) 


Points 1-3: 


m 


1-3 


f(p.u. +p,uj y.yj 1 

= (2jt) 6 | Pi v i yj (Xj-Xj)]. 


(4.161) 
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The particle density relationship at point 3 may now be determined by 
equating Eqs. (4.160) and (4.161) and solving for p 3> 

Equation (4.159) becomes: 

(Pj Uj +p 5 u 5 ) y t yj - 2p lVl y^ (Xj -x,.) = (p^ +p 3 u 3 ) ypr* - Zp^yJ (^ -x 3 ) (4.162) 

To solve for p 3 : 

combine terms of Eq. (4.162) to form coefficients of Pj.P 3 and P 5 
P l[ 2v l y l +P 5 U 5 y l y l = P 3 u 3 y l y l 


and then divide through by u 3 y^j 


P 5 U 5 1 Zp i V l (x 3' x 5 ) 

u. 


(4.163) 


4.4.3 Particle Limiting Streamline Point Solution 


Consider the characteristic net for a point on a particle limiting stream- 


line 
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Applying the principle of conservation of mass to species j in the above 
system, we may write 


• • « 

= mg _ 3 + « l 3_4 • (4.164) 

By definition of a particle limiting streamline, the entire mass of 
particle species j is contained in the streamtube formed by the particle 
limiting streamline for species j. Therefore ** 15.3 may be set equal to zero. 

Note that this condition applys only to particle species j; other particle species 
may pass through the particle limiting streamline for species j. Equation (3.164) 
then becomes 

m 3 4 = m 5 _ 4 . (4.165) 

Each term in the above equation may be solved for by applying Eq. (4.137) be- 
tween the proper points. 


Points 3-4; 


m 


3-4 


= (2 ir)° 


6 6 6 6 
^3 u 3 +f> 4 U 4 ) iy 3 y 3 ‘ yr 4 y 4 ) (P 3 V 3 +P 4 V 4 ) {y 3 + / 4 ) 


(x 3 - x 4 ) 


(4.166) 


Points 5-4; 


m 


5-4 


= ( 2 jr) 


[ 0 ) 5 U 


+P 4 U 4 ) ( y s Y s - y 4 y 4 ) (p 5 v 5 +p 4 v 4 ) (y c 


•y 4 ) 


( x 5 - x 4 ) 


(4.167) 


The particle density relationship at point 3 may new be determined by equating 
Eqs. (4.166) and (4.167) and solving for p 3> Equation (4.165) becomes: 


f f f f 

(P 3 U 3 +p 4 u 4 ) (y 3 y 3 - y 4 y 4 ) - (p 3 v 3 + p 4 v 4 ) (y 3 + y 4 ) (x 3 -x 4 ) = 

(p 5 « 5 +P 4 U 4 ) (>^5 - y 4 y 4 ) - (P 5 v 5 u P 4 V 4 ) (y| + y\) ( x 5 - x 4 ) • (4.168) 
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To solve for p^: 
expand Eq. (3.168) 

P 3 U 3 (y 3 /5 " y 4 V 4 ) ' P 3 V 3 (y 3 + y 4> (x 3 ' x 4 ) + P 4 U 4 (y 3 y 3 " *4? 4* 

~ P 4 V 4 (y 3 + y 4> (x 3 ' X 4 } = ^5 U 5 (y 5 y 5 ' y 4 y 4* * P 5 V 5 (y 5 +y 4 )(x 5 

6 6 6 6 

+ P 4 u 4 (y 5 y 5 - y 4 y 4 ) - P 4 v 4 (y 5 +y 4 ) (x & -x 4 ) ; 


combine terms to form coefficients of p^,P 4 and p^ 

(y 3 y 3 ' y 4 y 4> - v 3 (y 3 + y 4> (x 3 * x 4>] = ^["s (y 5 y 5 “ y 4 y 4> 

- v r (y 3 + y 4 ) ( x 5 - x 4 >] P 4 j u 4 ^ y 5 y 5 * y 3 y 3^ ' V 4^ y 5 +y 4 ) < x 5 " x 4 > 

- (y 3 + y 4 ) ( x 3 ■ x 4 )j| * 

and finally solving for p^ 

P3 * ——6 T —.6 - « : . : P5[ u 5 <y 5 y 5 - ' , 4' r 4» 

u 3 (y 3 y 3 - y 4 y 4> ‘ v 3 (y 3 + y 4* (x 3 * x 4* [ 

" v 5 (y 5 + y 4> (X 5 - x 4>] + " y 3 y 3> ’ v 4 [ (y 5 + y 4> < x 5 ' X 4 } 


- (y 3 + y 4 ) 


(x 3 - x 4>]| 


(4.169) 


4.4.4 Upper Boundary Point Solution 

The upper wall point and free boundary point are the two types of upper 
boundary points encountered in a typical gas-particle flow problem. Consider 
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the characteristic net for a point on an upper boundary 



Applying the principle of conservation of mass to species j in the above 
system, we may write 


”4-4 = ”4-3 * ”4-4 * (4.170) 

Unlike the assumption made in the lower boundary solution, the particle 
streamlines are not restricted to run parallel to the upper boundary in the 
vicinity of the boundary. Therefore, the particle mass flow rate m£ , is not 
necessarily equal to zero. Physically speaking, 3 may be considered to 
be equivalent to particles "sticking" to the upper wall or passing through the 
free boundary. 


The particle density relationship at the upper boundary is derived by 
following the same mathematical procedure used in the solution for the particle 
density relationship at an interior point. The final result is 


P 3 = 


A 

u 3 (y 5 y 5 ' y 4 y 4> * v 3 [ (y 5 + y l ) < x 5 + <4 +^3) ( x 3 - x 4 )l 

•> L -11 1 J 


- v r 


- v. 


(yg +y 4 ) ( x 5 - x 4 ' - (y^+y*) (x 5 -x 3 ) 


(y 5 + y 4 ) < x 5* x 4> * (y 4 +y 3 ) ( x 3" x 4 > 
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u 5 ,y 3 y 3 • 


, 6 6v 

4 ( y 5 y 5- y 3 y 3 ) 


(4.171) 
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5. SHOCK CAPTURING OPTION 


A shock capturing option has been included in the RAMP code to permit 
the computation of two-phase (gas-particle) nozzle flows with embedded 
shocks. The gas phase may be treated as either an ic^al or real gas with 
equilibrium/frozen chemistry. The particle phase is treated as a discrete 
distribution of particle sizes. The two phases interact only in the 
transfer of momentum and energy. There is no mass transfer between phases, 
i.e., the particles do not evaporate, and the gas does not solidify or 
condense. 

The inviscid steady state conservation equations for a gas-particle 
mixture are given as follows in vector notation: 

Gas 


Mass 

Momentum 

Energy 

Particles 

Mass 

Momentum 

Energy 


7»pq ” 0 

NP 

P*f(q) pq] +vv + ^2 P^ A-* A q^ =0 

j-1 

NP 

7-(pq H) - P^ A^ (B^ - q • Aq-*) “ 0 


7.p j q - 0 

V • [(q) J p - * q ■*] - p^ Aq ■* - o 

V'(p q h J ) + p J A J (B J - Aq . Aq J ) “ 0 


(5.1) 
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These equations, expressed in compact form in nozzle coordinates, become : 


+ -|? In 5 <F-«KP)l if* *H-0 (5.2») 

and, for axis solution: 

H + i 4 - <f " cP)+e <5 - 2b) 

where B, P and H are the following arrays of variables: 


Gas 
E - 



F * 


Pv 

puv 

pv 2 + 
pvH 



H 


Au^ 

Xo^A^ Av^ 
£p j A J (B j - 



Particles 


E 


pV 

P V 2 

pW 

P^u^ 


pV^ 

p^u^v^ 

p^v^ 2 

pW 


H 


0 

-p^A^ Au - ^ 

-p^ Av^ 

P J A J (B J - Au J2 - 



( 5 . 4 ) 


and where 

g 0 2D planer 

1 axisymmetric 

m 1 r or y momentum eq. 
0 all others 


The nozzle coordinates x and D » r/R or y/R are based on an expanding 
scale In the r or y direction proportional to the nozzle radius (or half- 
width) R. This permits a uniform distribution of grid points in the r or y 
direction. 
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Equation (5.2) is numerically integrated by the McCormack predictor/ 
corrector algorithm which is second-border accurate. This algorithm is given 
as follows: 


Predictor Step 


E ®E n +— (E“ - E R ) 

m m ^n An m+l m 


R \ 


Ax ,6 ,_n 

An m+l 


6 k Ci> 


(ni -6An) (*” - 6 <p, 

m ' m 1 


:»] 


- 6* 


Ax 


r° An 


(p“ . - P n ) - H n Ax 
m+l m' m 


(5.6a) 


and, on axis (m = 1) 


n- E "-7^[ (F "-' :p 2 ) - < p ; - «;>] -»i 


Ax 


(5.6b) 


Corrector Step 


C 1 ■ r |< E » + V + ^ % < E „ - Vi> 

[>4 


Ax 

6An 


+ Mr,) (F m - 6<P m ) - n; (Vi - «<p, 


»-i>] 


“ 6ic r— (P - P .) - H Ax 
An m m-1 m 

and, on axis (m *■ 1) 


(5.7a) 


_n+l 1 i n , - 2 Ax . - . - - , , - , 

E ■ 2 j E m + E » - -^TT ST l(P 2 ' kP 2 ) ' <P 1 ‘ KP 1>1 - H 1 ** 


(5.7b) 
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The nozzle wall boundary la treated using the HOC option In RAMP for 
the gas phase. The particle phase upper boundary does not necessarily 
follow the nozzle wall boundary, and, hence, must be tracked for each 
particle specie by integrating the following equations of motion: 


x • A Au 

r (or y> » AAv (5.8) 

These equations are Integrated with two steps along with Eqs. (5.6) and (5.7) to 
find the limiting streamline for each particle specie: 

Predictor Step 



R° + A n 

P W 



Ax 


2 < u »> 



Corrector Step 


(5.9a) 



£ | R n + R + A Av 
2 ( p p mm 


Ax' 


m 


2 ( u _) 

m 


Ax 


u 


m 


(5.9b) 


The flowfleld properties are obtained by decoding the computed flux 
variables E obtained by integrating Eqs. (5.6) and (5.7). The decode procedure 
for the gas phase depends on the type of chemistry involved. For both ideal 
and real gas decodes the radial (vertical) velocity v and total enthalpy H 
are obtained by 


v 


n+1 


e 3 /e 


1 


H 1 * 1 
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In the Ideal gas decode, <he axial velocity u is obtained from: 



( 5 . 10 ) 


The real gas decode procedure makes use of real gas chemistry tables in 
RAMP, These tables, by a look-up procedure, recover temperature, pressure, 
density, Mach number, gas constant and ratio of specific heats from an input 
of total enthalpy, entropy and velocity. The shock capture option real gas 
decode uses Newton iteration to solve the following equation: 



(5.11) 


where the gas constant, R, and temperature, I, are obtained from the real 
gas tables from known values of entropy, S, total velocity, V, and total 
enthalpy, H. Following convergence of the iteration on velocity, V, the 
axial velocity, u, is obtained from: 


n+1 

u 


a 



(5.12) 


For both Ideal and real gas decodes, the density, p, and pressure, p, 
are then obtained from: 


n+1 

P 


Ej/u 


n+1 


(5.13) 


n+1 

P 


n+1 

u 


(5.14) 


5-5 


LOCKHEED- HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC TR D867400-I 


The particle phase decode Is: 


ttf 1 

(U J ) 

- e^/eJ 

(5.15) 

. n+1 

. 2 . 


<p J > 

- <eJ) 

(5.16) 

. n+1 

(v J ) 

- eJ/eJ 

(5.17) 

. n+1 

(h J ) 

- rJ/eJ 

(5.18) 


The Shock Capture option was checked out by comparison with the 
existing streamline/ normal method-of-characteristlcs option. Three test 
cases were considered: (1) SSME nozzle flow, single phase, ideal gas, (2) 

SSME nozzle flow, single phase, equilibrium chemistry, and (3) two-phase 
ideal gas with solid particles. Computed Mach number and pressure distribu- 
tions are shown in Figs. 5-1 through 5-3 for case 1 (SSME nozzle, ideal 
gas). Axial distributions are shown in Figs. 5-1 and 5-2, and exit radial 
distributions are shown in Fig. 5-3. Generally good agreement is shown 
between the Shock Capture and Streamline/Normal options. Similar results 
are shown in Figs. 5-4 through 5-6 for case 2 (SSME nozzle, equilibrium 
chemistry). Computed Mach number and pressure distributions along the exit 
.adlus are shown in Fig. 5-7 for case 3 (Space Shuttle SRM, two-phase, ideal 
gas). Solid particle densities, velocities and temperature exit radial 
distributions are shown in Figs. 5-8 through 5-10. In all cases, good 
agreement is shown between the Shock Capture and Streamline/Normal options. 

On the basis of the demonstrated good agreement with the Streamline/ 
Normal option, we conclude that the Shock Capture option is operational and 
error free. 
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Fig. 5-4 Mach Number Distribution Along SSME Nozzle Axi 
Equilibrium Chemistry Solution 






Pig. 5-6 Mach Number and Pressure Distributions at SSME Nozzle Exit 
Equilibrium Chemistry Solution 
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Fig. 5-7 Mach dumber and Pressure Distributions Along Nozzle Exit Radius, 
Two Phase/Ideal Gas Case 



Algori cri.j Option 



Fig, 5-8 Solid Particle Density Distributions Along Nozzle Exit Radius 
Two Phase/Ideal Gas Case 




Algorithm Option 



Fig. 5-9 Solid Particle Velocity Distributions Along Nozzle Exit Radius 
Two Phase/Ideal Gas Case 





Fig. 5-10 Solid Particle Temperature Distributions Along Nozzle Exit Radius 
Two Phase/Ideal Gas Case 
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6. EXPANSION CORNER — PRANDTL-MEYER FAN 

In some cases the flow may be required to negotiate a sharp expansion 
turn. The problem becomes two dimensional at a sharp corner (it is impos- 
sible to conceive of an expansion corner on an axis of symmetry) and may be 
treated with a Prandtl -Meyer expansion. 


q+dq 



Fig. 6-1 - Nomenclature for Mach Wave Analysis 


Since a Mach w?ve will support pressure changes only in a direction normal 
to itself (Ref. 4iv 


dv = qd G 
du = dq 


du 

dy 


tan i u = - — ■ ■ 

^ M 2 - 1 


6 - 1 


LOCKHFED-HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC TR D867400-I 


or 



( 6 . 1 ) 


The solution to Eq. (6.1) is a straightforward numerical integration for 
the case of a known final velocity (free -boundary case). If the turning angle 
is known, however, and the final velocity is not known, an iterative solution is 
necessary to d etermine the upper limit. 


% 



f - * « 1 2 > * 0 


( 6 . 2 ) 


In the mesh construction to be discussed later a fan of rays must be gen- 
erated to allow a numerical description through a large turning angle. The turn- 
ing angle is subdivided into a number of small turns, each of which is integrated 
numerically. Corresponding to each of these small turns is a Mach wave or 
characteristic line emanating from the corner. 
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7. NUMERICAL SOLUTION (MESH PCINT CONSTRUCTION) 

The calculations described previously are point or small region solu- 
tions, This section presents the mesh construction and calculation procedure 
required to describe the entire flow field in a typical gas-particle flow problem. 
The general principles adopted in the numerical solution are similar to those 
used in the method-of-characteristics solution*?. The major difference is in 
the technique of mesh construction; here, the calculation proceeds along 
normals to the streamlines instead of along characteristic lines. 

To begin the problem of describing the entire flow field all necessary 
boundary conditions must, of course, be supplied. In addition, a start line 
whose properties are completely defined must be designated. 

Figure 7-1 lilustrates a flow f ; eld in which there are no discontinuities 
and in which the mesh construction it terminated when the region of interest 
has been computed. 

Figure 7- 2a pre&-^n*:s the mesh construction required to solve for the 
properties of any interior point in Fig. 7-1. The J ^ne is considered to be 
the known normal line. The J-line may oe the input startline or the line just 
calculated. The K-line is the new norma, line to be calculated. Calculation 
always starts from the lowest point and moves upward along the K-line. The 
first point on the new normal is located by extending the right- running char- 
acteristic from a properly selected point (based on desired mesh size) on 
the preceding knovm normal to intersect the lower boundary based on the 
flow properties of that point. The rest of the points including the ^ast point 
(boundary print) of the normal are located by extending the normal to the 
local streamline to intersect the naxt streamline (A boundary is also a 
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streamline.) Once the location of the new point is known (point 3 in Fig. 7- 2a), 
the characteristic line(s) and particle streamline(s) can be made to intersect 
the preceding known normal (points 1, 2 and 8 in Fig. 7- 2a). The properties 
of points 1, 2 and 8 on the known normal can be interpolated from the known 
points. The properties at point 3 can then be calculated by solving the applic- 
able compatibility equations. For an interior point, both the right-running 
and left-running characteristics are used to compute the velocity and the flow 
angle at the new point (point 3 in Fig. 7- 2a). For a boundary point, either 
the flow angle or the velocity can be found with the boundary condition; there- 
fore, only one physical characteristic is needed depending on the relative 
location of the boundary with the main flow. For shock points and Prandtl- 
Meyer expansion calculations, the mesh construction is appropriately modified 
to accommodate each individual case. The general method of computation, 
however, remains the same as described above. 

The mesh construction and calculation procedure for an interior point 
solution forms the basis for the solution of the flow properties at all other 
types of points in the flow field. For that reason, the mesh construction 
and calculation procedure for an interior point will be first described in 
detail, then, to further aid in the understanding of the numerical solution, 
presented in flow chart form. Using the knowledge gained in the discussion 
of the interior point solution, the calculation procedure for each of the re- 
maining cases will be described. Any important problems that may arise, 
or other points of interest involved in the calculational procedure will be 
discussed. 

7.1 INTERIOR POINT 

This section presents the details of the numerical solution for the gas 
and particle flow properties at an interior point of the flow field being analyzed. 
The solution is outlined in a step-by-step procedure and is for the case when 
particles are present at all points which enter into the calculation of the flow 
properties at the new point. 
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Consider again Fig. 7-2a and assume that the flow properties are known 
along the J-line. The J-line can be the input start line or the line just calcu- 
lated. The point 7K will also be a known data point since it is the last point 
calculated on the K-line. 

Before the calculation procedure can be started a decision must be 
made to use either the enthalpy-entropy-velocity form or the pressure- 
density-velocity form of the compatibility equations to solve for the proper- 
ties at the new point. This decision is based on the type of chemistry being 
considered in the flow field. If finite rate chemistry is being considered, 
the pressure -density-velocity form of the compatibility equations is the more 
convenient form to be used. If equilibrium or frozen chemistry is being con- 
sidered, the enthalpy-entropy-velocity form of the compatibility equations is 
the more convenient form to be used. In either case, the method of calcula- 
tion is the same. For purposes of discussion, it will be assumed that the 
chemistry of the flow field being analyzed may be considered to be in chem- 
ical equilibrium. 

Step 1: For the first pass of the solution, the properties at 3K are set 
equal to those of the base point 5J. 

Step 2: The physical location (x^, of the point 3K is found by inter- 
secting the normal from the point 7K with the projection of the streamline 
from the point 5J. 

Step 3: The intersection of the characteristic lines with the J-line is 
now found. These lines will in general fall between two known data points 
and an interpolation scheme is employed to find the flow properties at these 
points. The interpolation is of the following form: 


P 1,2 = P 5> 4 + F < P 6,5- P 5,4> 
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where 


F = interpolation factor. 

For the gas flow properties P assumes the following 


with 


P =q,9,S,H 


T =f(q,S,H) 
p * f (q, S, H) 


v - f (T,H) 
C p s f (T,H) 

Pr = f (T, H) 


tabulated 


while for the particle properties P assumes 

P = J,v*,h^,p^ 


with 


^ = f(h^)j ^ tabulated 


Re 


1,2 


2r l,2 p l,2 M l, 2 

v \, 2 


fj 2 = 2 ^bulated 

C l,2 = f « Re i,2- Pr l,2> Ret - 12 
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(3.46) 





Pr 


1.2 


(3.83) 


and. 


B, 


1.2 


P 1 . 2 


R 


- 1 


q l,2* Aq l,2 " q l,2 


A< !i,2 + ! C i,2< T i,2- T l,2» 


1.2 


+ 


3a 


1.2 


A i,2 in i,2 r i f 2 





(3.86c) 


Note that the points 1, 2 and 8 are not necessarily as shown in Fig. 7- 2a. 
Their locations can vary along the J-line depending on the mesh size and the 
Mach number in the vicinity of point 3. Figure 7-2b presents one example 
when the point m + 2 is an upper boundary point and the point 1 ! is outside the 
flow field under consideration. Instead of the point 1* the point 1 is used in 
the calculation. However, the properties at point 1 are not available nor can 
they be interpolated because point n + 2 is not yet known. For t i_.*s particular 
case, the flow properties at point 1 are assumed to be identical to those of 
point m + 2 so that the properties at point 3 can be approximated. The same 
method is used to calculate point nfl. The boundary point n {- 2 of the In- 
line can then be calculated without problem. Once the boundary point is deter- 
mined the calculation goes back to point 1 and the right running characteristic 
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Fig. 7-2b - Mesh Construction for an Interior Point 


intersection is made using the two boiuu.i r points. Point 1 is then solved 
for and the solution proceeds to the plume boundary. 

Step 4: The particle density at point 3K is calculated using the relation- 
ship: 
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Step 5; The particle streamline intersections with the J data line are 
now located (point 8 of Fig. 7- 2a) ucing the particle streamline relationship: 


yj - J 
_3 ^8 

x 3 “ x 8 


(4.42) 


at 


Step 6: The gas and particle properties at point 8J are now found by 
employing an interpolation scheme as described in Step 3. 

S tep 7: The coefficients for the particle compatibility relations are 
now calculated from Eq. (4- 150a). 


C zi 


j _ ( 8, 3 ~ 8 , 3 ) j 

8,3 " j A 8 

u 8, 3 


‘8, 3 


C 3 


8,3 


_ ( v 8,3 - J 


u' 


‘8,3 


8,3 


C4 


3 A 8,3 8,3 Vi 8,3 " "8,3 


(Ti - - T q , ) + 


3o 


8,3 


8,3 


m 8,3 r 8,3 


e 8,3 u 8,3' “8,3 A 8, 3 


U' 


J 

8,3 


Step 8: The new particle properties at point 3K are then found by ex 
panding the difference Eqs. (4.148), (4.149) and (4.150) to obtain 

u 3 = u 8 + C2 8,3 Ax 8-3 
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H = ^ + ^8,3 Ax 8-3 

h 3 = h 8 + C 4 8, 3 Ax 8-3 
with the remaining properties defined as 


= f (h)^ tabulated 


Re ) Zf 3 P 3 S 

3 - K, 


fj = f(Re^) tabulated 


= f(Re^, Pr 3 ) Ref. 12 


&J 

A 3 


9 f 3^3 

1 


m 1 


J 

3 '*3 


TT 


(H) 


and. 


g! C 

cl . JJL 

f 3 Pr 3 


B l, * c /r\ - 1 +f c 3< t 3- t : 

3 P 3 ' 3 1 


H ,T i , 4 -“i T t]| 
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Step 9: The coefficients for the gas streamline compatibility relations 
are now calculated 


C 2 


P 5,3 


- 1 


R 


5,3 


NP 


P 5,3 A 5,3 B i, 


AL 


5-3 


5 ’ 3 J P 5,3 T 5,3 q 5,3 


and 


NP 


5.3 j=l 


P 53 A 5,3 < u 5,3 " U 5,3 ) cos6 5,3 + (v 5 ,3~ v 5,3* Sul6 5,3 1 AL 5-3 


Step 10: The rate of change of entropy in the direction of the gas stream- 
line and the entropy at the point 3K are calculated in the following manner: 


AS 


5-3 



(4.143) 


and the entropy at point 3K is 


S 3 = S 5 


Step 1 1 : The rate of change of enthalpy in the direction of the gas stream- 
line and the enthalpy at the point 3K are calculated in the following manner: 


AH 5-3 “ T 5,3 AS 5-3 



(4.142) 
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and the enthalpy at point 3K is 

H 3 = h 5 + T 5>3 AS 5 3 - . 

Step 12: The remaining gas properties to be calculated at point 3K are 
the velocity and flow angle. Expanding equation (4.146) in terms of both the 
right and left-running characteristics and solving for the velocity gives 


e i- 0 2 +Q l < ll +Q 2 < i2- B r B 2 + CI l< H 3“ H l ) + C ¥ H 3‘ H 2> + G l +G 2- C1 l- C1 2 

q 3 = ^ 


1 + *2 


The new flow angle is then 


0 3 = 6 2 + Q 2 q 3 - Q 2 q 2 + B 2 " C1 2 ^ H 3 ' H 2^ ” G 2 + C *2 


Step 13: The properties T 3 , Py Vy Cp , Pr 3 , A^, and are recal- 
culated to reflect the new gas properties. The solution is then checked to 
determine if the desired accuracy has been obtained. If the solution has not 
converged, return to step 2 and repeat the process. The flow of the numerical 
solution is summarized in Chart 7- 1. 
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Chart 7- 1 - Flow Chart of the Calculation Procedure for an Interior Point Solution 
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7 . 2 LOWER WALL POINT 

Figure 7-3 presents the mesh construction required to solve for gas and 
particle properties at any lower wall point in Figure • * 1. 


Normal 


Particle 

Streamlines 



Right-Running 
Characteristic Line 


VM7mr 

3K 

Gas Streamline 
Along Lower Wall 


Fig. 7-3 - Mesh Construction for a Lower Wall Point 


The calculations! procedure followed to determine the flow properties 
at a lower all point is very similar to that of an interior point solution. The 
lower wall point is always the first point calculated on a K-line. Before the 
point 3K can be located, jpoint 1J must first be selected as a function of the 
desired mesh size. Once point 1J is located and its properties determined 
using the interpo? ation scheme employed for an interior point; the physical 
location (x^, y^) of point 3K is found by intersecting the right-running charac- 
teristic from point 1 J and the projection of the streamline from the point 5J. 
The particle density at point 3K is then calculated using the relationship 



P J 5 u J 5 - 2 pj vJ (x 3 -x 5 ) 




(4.163) 


and is based on the assumptions that: (1) the lower boundary runs parallel to 
the X-axis; and (2) the particle streamlines run parallel to the lower boundary 
when in the vicinity of the lower boundary. 
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Finally, since the flow angle along the lower wall is known, only the ri^ht 
running characteristic line is needed to calculate the velocity at the lower wall 
point. Solving for velocity in Eq. (4.146) gives 


[ 0 l' 6 3 +Q l q l * b 1 + CI 1 (H 3 -H 1 ) + Gj-ClJ 



Once the velocity is found, the calculational procedure is continued as 
in the case of an interior point 

7.3 UPPER WALL POINT 

Figure 7-4 presents the mesh construction required to solve for the gas' 
and particle properties at any upper wall point in Fig. 7-1. 



Fig. 7-4 - Mesh Construction for an Upper Wall Point 

The calculational procedure followed to determine the flow properties at 
an upper wall point is almost identical to that of an interior point solution. The 
exception is that the flow angle on the upper wall is known andjti.erefore, only 

the left-running characteristic line is needed to calculate the velocity at the 
upper boundary point. Solving for velocity ir Eq. (4.146) givest 
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® 3"®2 + ^ 2^2 ’ B 2 + ^* 2 ^ H 3 _H 2 ^ + ^’ 2”^*2 

q 3 = Q, 


7.4 FREE BOUNDARY POINT 


Figure 6-5 presents the mesh construction required to solve for the gas 
and particle properties at a free boundary point. Note that the mesh cons .auc- 
tion is the same as that for an upper wall point. 



F ree Boundary 


Left-Running 
Characteristic Line 


Fig. 7-5 - Mesh Construction for a Free Boundary Point 

The calculational procedure followed to determine the flow properties at 
a free boundary point is like that of an upper wall point solution. The one dif- 
ference is chat the velocity and not the flow angle at the free boundary point is 
known. Applying Eq. (4.146) to the left- running characteristic and solving for 
the flow angle at the free boundary point gives 

0 3 = e 2 + Q 2 <q 3 -q 2 )*B 2 -CI 2 <H 3 -H 2 )-G 2 + Cl 2 . 

The velocity at the free boundary, though assumed to be known, is not 
easily obtained. The velocity at the free boundary must be solved for itera- 
tively using the relation 
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An interative solution is required because the gas properties R, T and P 
are functions of a variable y. If an ideal gas was being considered, y would be 
constant and may be solved for directly, 

7.5 PARTICLE LIMITING STREAMLINE POINT 

Figure 7-6 presents the mesh construction required to solve for the gas 
and particle properties at any particle limiting streamline point in Figure 7-1. 



Fig. 7- 6 - Mesh Construction for a Particle Limiting Streamline Point 

The calculational procedure followed to determine the flow properties at 
a particle limiting streamline point is somewhat diLerent from th't of an inter- 
ior point solution. The differences in the calculational procedure are brought 
about by the fact that the J and K normals are constructed between a gas stream- 
line and a particle limiting streamline rather than between two gas stream- 
lines. 


The physical location (x^, y^l of point 3K is found by intersecting the 
normal from the point 7K with the projection of the particle limiting stream- 
line of species J from the point 5J. Once point 3K is located, the gas streamline 
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intersection with the J data line must be located (point 9 of Figure 7-6) using 
the gas streamline relationship 





(^• 41 ) 


The properties at this point are found by interpolation between the known points 
4J and 5J. 


The particle density at point 3K for all particle species except particle 
species J is determined in the same manner as for an interior point. The par- 
ticle density of species J is calculated using the relationship 


*3 = 


1 


u 3 (y 3 y 3 ‘ y 2 y 2 } ' v 3 (y 3 + y 2 )(x 3' X 2 ) 


h 


P 5 | U 5 (y 5 y 5 -* 2 * 2 * ' 


(x 5 - x 2 ) 


] +P 2 J u 2 (y 5 y 5 • y 3 y 3> ‘ v 2 [ (y 5 + y 2> (x 5 ' x 2 ) ' (y 3 + y 2> (x 3 * 


(4. 


and is based on the fact that the entire mass of particle species J>is contained 
in the streamtube formed by the particle limiting streamline for species J. 


Once the particle densities nave been found, the calculation procedure is 
continued as in the case of an interior point. 
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7.6 PARTICLE LIMITING STREAMLINE-BOUNDARY INTERSECTION 

Figure 7-7 presents the mesh construction required to solve for the gas 
and particle properties at the particle limiting streamline -boundary intersection. 



Fig. 7-7 - Mesh Construction for a Particle Limiting 
Streamline -Boundary Intersection 

Assume that the particle limiting streamline point (n, K) has just been 
completed and the next point to be calculated is that of a free boundary point 
(n + 1, K). Under these circumstances, the free boundary point is calculated 
in the normal manner, and then upon completion of the free boundary point a 
comparison is made of the physical location of the two points. Ify n+ j< y^, 

an intersection of the particle limiting streamline and free boundary is indi- 
cated. The physical location of (n+2, K) is four, * intersecting the projec- 
tion of the streamline from the point (m + 1, J) am ..he particle limiting stream- 
line between the points (m, J) and (n, K). The gas and particle properties at the 
point (n+2, K) are then found by interpolation between the points (m, J) and (n, K). 

Once the properties at the point (n+2, K) are found, the pomt (n+ 1, K) is 
assumed to no lo-^^r exist. The point (n+2, K) is used as the base point the 
next time an upper boundary point is calculated. 

The description of the entire flow field now proceeds onto the calculation 
of the lower boundary point of the new K-line. 
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7.7 EXPANSION CORNER- PRANDTL- MEYER FAN 

Expansion corners are found when the slope of the solid boundary has a 
discontinuity in such a way that the flow must negotiate a sharp turn away from 
the oncoming flow or when the flow issuing into an open space is under- expan- 
ded in a nozzle or a flow channel. The expansion takes place either in a chan- 
nel or at the channel exit. The amount of turn the flow must make and the fi- 
nal flow velocity can be found with the aid of the Prandtl-Meyer expansion re- 
lation, Eq.(6.1),and the boundary conditions. 

After the total turning angle is found, it is subdivided into a number of 
small turning angles. The flow properties at the expansion corner correspon- 
ding to each of the small turns are calculated and the angle' s corresponding 
characteristic lines emanating from the corner are constructed. The final 
line of the expansion is a streamline rather than a characteristic line and, 
normally, the Mach angle corresponding to this streamline is substantially 
greater than the small turning angle. An extra point is inserted to fill the 
gap after the expansion corner computation has been completed. Note that 
this final expansion line, a streamline, can be a solid boundary or a free 
boundary. 

In the example shown in Figure 7- 8, the turning angle is subdivided into 
four small turns. The corresponding characteristic lines emanating from the 
corner are designated as A, B, C an.; D in Figure 7-8. The final expansion 
line E is a streamline corresponding to the last corner point (d, J). The pro- 
perties at the corner corresponding to every small turn are stored in (a, J), 

(b, J), (c, J) and (d, J) arrays. The properties at the corner before the ex- 
pansion takes place can be found as usual and are stored in (m, J) array; its 
characteristic line is designated as M in Figure 7-8. 

A weak shock is initialized at point (n+4, K) for the upper corner expan- 
sion point or (n-4,K) for the lower corner expansion, which is done simply by 
defining a double point at those points with identical flow properties. The initial 
shock angle is assumed equal to the Mach angle. 
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Fig. 7- 8a* Mesh Construction for an Upper Expansion Corner 
Point 



Fig. 7- 8b- Mesh Construction for a Lower Expansion Corner 
Point 
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7.7.1 Upper Corner Expansion 

As shown in Figure 7- 8(a), the point (n-1, K) is considered to be com* 
plete. A left-running characteristic is drawn from the point (n-l.K)to intersect 
the right- running characteristic from point (m, J) at M, which is an interior 
point. The properties can be calculated in the usual manner. A normal is drawn 
from point (n-1, K) to intersect the right- running characteristic from point (m, J) 
at n. The properties at this point are then interpolated between points (m, J) and 
M. Similarly, the properties at points (n+1, K) through (n+4, K) can be calcu- 
lated. Point (n+5, K) is calculated by the same technique provided that the point 
E is treated as an upper boundary point. If the distance between points (n+4, K) 
and (n+ 5, K) is greater than other meshes in the expansion fan, one extra point 
(n+6, K) is inserted and its properties are found by interpolation between points 
(n+4, K) and (n+5, K). 

7.7.2 Lower Corner Expansion 

This case is shown in Figure 7- 8 (b). When an expansion at the lower 
boundary is detected, the normal line starting from the corner point (m, J) is 
completed first. The flow properties at the corner and the corresponding 
characteristic line for each small turn can be found by the same method de- 
scribed in the preceding paragraph. Point (n, K) is the intersection of the 
right- running characteristic from point (m+1, J) and the left- running charac- 
teristic from point (m, J); its properties can be calculated without problem. 

The flow properties at point A can also be calculated. A normal from point 
(n, K) and the left- running characteristic from point (a, J) intersects at (n-1, K); 
flow properties at this point are then interpolated between points (a, J) and A. 
Points (n-2, K) through (n-4, K) can be calculated by the same method. Point 
(n-5, K) is a lower boundary point which can be calculated with the technique 
used to find the point (n+5, K) in Figure 7- 8(a). An extra point may also be 
inserted between points (n-4, K) and (n-5, K) if needed. 
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8. NUMERICAL INTEGRATION OF THE CONSERVATION EQUATIONS 

Examination of the mass flow rate, momentum flux and energy flux 
through a nozzle and exhaust plume will yield considerable information about 
the nozzle performance and the status of the flowfield numerical solution. 
The momentum flux is used to determine the thrust produced by the nozzle 
while the thrust ratioed to the mass flow rate is a measure of the nozzle 
performance. Deviations of the fluxes in the conservative parameters from 
surface to surface (Fig. 8-1) are indicators of the validity of a numerical 
solution such as the one used in the gas -particle code. 



Fig. 8-1 - Schematic of Conservation Parameters in a Nozzle Plume 
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In terms of the coordinate directions, the incremental area is 

7 7 i/2 

dA = (2jrr)® |(dr) 4 + (dx) f 


( 8 . 2 ) 


The flow area in the direction normal to the surface, dL, is now 

dA = (2Tr) 6 {(dr) 2 + (dx) 2 }^ 2 e n 


(8.3) 


where the unit normal vector, © n » is defined by the geometry of Fig. 8- 2c as 


A a 

e n = cosip i + simp j ; 


(8.4) 


and the angular location of the unit normal vector relative to the x-axis is 

<p = - tan' 1 (8.5) 

The quantity of mass flowing through the incremental area is found from 
the relation 


dm = p|^|dAe q *e n ; 


( 8 . 6 ) 


where the unit velocity vector, 0 is defined by the geometry of Fig. 8-2b. 
The relation for m is now 

m = p | q JdA {cos0 i + sinft • (cos <p i + sirup j } . 


Performing the indicated dot product the relation reduces to 


m = p |q|dA {cosd cos<^> + sin0 sin<p| . 


The trigometric identity for sine-cosine function is 


cos0 cos(p + sin0 simp = cos (0-<0 . 
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Substituting, the incremental mass flow is 

dm = p J q |dA cos(9-<p) . 

The total mass flow through the surface is simply 

m = ^ dm 


(8.7) 


( 8 . 8 ) 


for all the incremental flow areas. 

So far the discussion has been directed toward a general expression 
for the mass flow rate. For gas -particle flows the surface geometry is the 
same. The only change is in the definition of the flow variables. The ex- 
pression for the respective phases are as follows: 


Gas Phase 


t.ll 

Particle (j particle) 


m = p|q dA cos(9 -<?) ; 


= p^ Jq^ldA cos(0^ -cp) ; 


(8.9a) 


(8.9b) 


Mixture 


NP 

“m = m + ^m-* 

j=l 


(8.9c) 


8.2 MOMENTUM 


As with the mass flow relations, a general expression for the momentum 
flux will be developed. This expression will then be extended to include the 
treatment of gas-particle flows. 
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The general expression for the momentum flux is 


= jp“ + P q qj • dA 


( 8 . 10 ) 


where 6 is the Kronecker delta and dA is given by Eq. (8.3). Substituting for 
dA and performing the indicated algebraic operations, Eq. (7.10) becomes 


= JPn + p q q .nj 


dA 


( 8 . 11 ) 


where dA is now defined by Eq. (8.2). Equation (8.11) yields an expression 
for each coordinate direction in the equation in terms of the scalar compo- 
nents given by the expressions 


= |p(cos<pi + si n<pj) + p(ui + vi) (u? + v j) • (cosipi + sin</?j)j. dA ; 


- [ 


M 




A A AA 

• i ii a \ 


A ) 


M = |P(cos</?i + sinCPj) + p(u ij + uvij +uvji+v ji) . (cosipi + sinx^j) j dA; 

-a 1 A, A 2 A A A 2 A ) 

= lP(cos<pi + simpj) + p(u cosxpi + uv cos</>j + uvsinc pi + v sin<^j)j dA 

The scalar equations in the coordinate directions are then 

2 S _,J 


M 


M 


and 


M 


M 


I 2 I 

= jP cos (p + p { u cos (p + uv si tup) j dA , 

( 2 ) 

= |P sinV? + p(uv cos <p + v sin</?)j dA; 


Substituting for the velocity components 


M 


M 


= Jp cos ip + pq^ (COS0 cosO cos (p + cosO sinO sint/?)j dA 


and 


( 2 ) 

M = ,P si tup + pq (cosO sinO cos + sin0 sin0 sin<^)| dA . 
M f I ) 
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Factoring common terms and substituting the trigometric ider._ity as before 
the final result for the incremental momentum’ flux is 


|p cosiP 4 /)q^ cos(0 - <p) cos0 1 dA , (8.12a) 

|Psin<p + pq^ cos(9 -*/>) sinoj dA (8.12b) 

The only difference between the expression for the gas and particle 
contributions to the momentum flux is thct the particles do not contribute to 
the pressure of the mixture. The respective expressions for the momentum 
flux through the surface are: 


and 


M 


M 




Gas 


= jP cosip v pq*~ cos(9 - ip ) cos0 \ A 


2 \ 

KLj = P sirup + pq cos(0 -</?) sin0j A 

I 


Particles 0 *^ Particle) 


= p^ (q^) coe(0 -<P) cos(0^) .. 


= p^ (q^) cos(0 -<p) sin(G^) A 


M. 


Mixture 


M 


M 


m 

x 


M 


M 


m 


NP 

• “m + X>m 

X j =1 x 
NP . 

* m m + S m m 

r j=l r 


(8.13a) 

(8.13b) 


(8.13c) 

(8.13d) 


(8.14a) 


(8.14b) 
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8.3 ENERGY 

The general relation is 


E = (h + y q^) (8.15) 

Applied to gas and particulate phases the respective relations for a given 
surface are: 


Gas 

E = (h + \ q 2 ) rii 


Particles (j** 1 Particle) 


E^ = (h^ + 4 (q 2 ) ) m J 


Mixture 


NP 

E = E + Y] E j 
m 

j-1 


8.4 NOZZLE PERFORMANCE 


(8.16a) 


(8.16b) 


(8.16c) 


Performance parameters considered are the thrust produced and the 
associated specific impulse, I . The force vector is calculated a. the 
mixture momentum on the initial data surface (Fig. 8-1) plus the pressure 
force on the nozzle wall. Mathematically this is expressed as 

f T * “m, + / P w < 817 » 


The integral term in Eq. (8.17) is obtained I / numerically integrating the 
pressure force along the nozzle wall. Consider the elemental surface, uL, 
bounded by the nozzle wall points a and b on any two successive data 
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Fig. 8-3 - Schematic of Pressure Force Acting on the Nozzle Wall 


surfaces. The slop *e surface, dL, is given by 


= tan_1 • 


(8.18) 


However, the fores due to the pressure loading acts normal to the sur < 
wh : ch is given by the angular relation 90 + ip. The local surface area over 
which the pressure is acting is 


dA = (2ir r) 


I (dr) 2 + (dx) 2 


I 


|l/2 

\ 


(8.19) 


and the unit normal to the surface is 


_k a a 

= - si.Kpi + cos<pj , (8.20) 
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so that 


dA = dA n 


dA = (2ff)® j(dr) 2 + (dx) 2 j | - sintpi + cos(<p)j| 


(P.2!) 


Equation (8.21) can be generalized to reflect the loading increment to either 
an »oper or lower wall in the following manner. Define a parameter X which 
has a value of Jh 1 for an upper or lower wall, respectively. The orientation 
of the surface unit normal with respect to the x-axis is 


9>1 = (p^ + X 90 

From the sine and cosine trigometric identities, Eq. (8.22) reduces to 


( 8 . 22 ) 


a^d 


co8<jPj = - sin <p 


si n^>j = cos q> 


so that the unit normal vector is now given by 


n = cos<Pji + sin<^|j ; 


and the surface area term is now 


dA 


w 


= (2»r) 6 |(dr) 2 + (dx) 2 j^ 2 jcos(^)i + sinft^Jj! . 


The resultant force relation for the nozzle is then 

NS 

— X 

F. 


s +2[P(2#r) 6 j(dr) 2 + (dx) 2 J 1/Z 1008(^)1 + sin^jj] (8.23) 

1 i=l ' ' 


and the nozzle thrust is just F . The specific impulse, I , is then 

x sp 
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*.p * VSo, • < 8 ' 24 > 

8.5 PERCENT CHANGES IN THE CONSEP’ ' HON QUANTITIES 

Okie check on the validity of the numerical solution is to monitor the 
change in the conservation quantities mass flow, momentum and energy. If 
these quantitie ao net change relative to the initial surface values, then the 
mass, momentum and energy have been conserved; this is a good indication 
that the numeri -<tl solution is proceeding satisfactorily. However, since a 
numerical solutirn is employed some changes in the conservation quantities 
is to be expected. The objective then is to keep the error as small as possible. 
When large errors occur, the Input data, wall equations, mesh construction, 
etc., should be checked. 

The changes in the conservation quantities are found by comparing ‘ie 
integrated value at a given data surface relative to the initial data surface. 
These are as follows: 


Mass Flow Rate (Mixture): 


> • * 


%Am = 


(mjq - nij) . 100 


m. 


Momentum (Mixture) 


%a|m| = 


<1*^1 - l»*i 1 • 100 


Energy (Mixture) 


%AE = 


(E n - Ej) • 100 


(8.25) 


(8.26) 


(8.27) 


The individual changes for the respective phases are computed in a similar 
fashion as specified by Eqs. (8.25) through (8.27). One important item to 
remember is that when gas particle flews are being considered the 
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solution is couple 1 in the sense that momentum and energy is exchanged 
between the phases. This means that relatively large changes in the momentum 
and energy quantities for each phase can occur. What should be of concern 
for gas-particle flows is the change in L’.ie mixture conservation quantities 
(momentum and energy). The respective mass flows should be conserved 
since the assumption has been made that the phases do not exchange mass 
via chemical reactions, condensation, etc. 

Finally, the specific impulse for a nozzle can be calculated in one of 
several ways. One way has already been outlined where the thrust is computed 
as the sum of the initial surface momentum and the pressure force acting on 
the nozzle wall. The specific impulse is then calculated by dividing the thrust 
by the mass flow rate through the initial surface. Using the initial surface 
mass flow rate is a matter of choice and can be argued from either a pro or 
con view. However, use of the initial surface mass flow rate eliminates the 
use of one variable which would be numerically calculated. Computation of 
the specific impulse by the second procedure is as follows. 


The thrust generated by the nozzle is first obtained by calculating the 

momentum in the r >zzle axial direction (M^j , Eq. (8.14a)). The specific 

m 

x 

impulse is then calculated by dividing by the mixture mass flow rate. This 

gives a means to compare the I calculation previously discussed. The 

sp 

change in I is 

Sp (M k . /m - I ) * 100 

M ' m sp 
m 

%AI = * j ( 8 . 28 ) 

p sp 


where the term I is obtained 
sp 


from Eq. (8.24) . 
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9. CONCLUSIONS AND RECOMMENDATIONS 

A computer code has been developed which can be used to model the 
dominant phenomena which affect the prediction of liquid and solid rocket 
nozzle and orbital (as well as low altitude) plume flow fields. With a 
single computer run it is possible for the designer to predict a high 
altitude plume starting from the combustion chamber. This code is intended 
to be a user-oriented design tool that will allow rocket nozzle/plume 
calculations that can range from the simplest preliminary design 
calculations to final design predictions. 

The emphasis in developing this tool was to require very little user 
interface and data manipulation which has previously been required to 
adequately treat such influences as nozzle boundary layer, two-phase and 
variable 0/F transonic startlines, boundary layer effects on particulates 
entering the boundary layer, and the transition from continuum to free 
molecular flow. Development of the program to this stage also allows the 
RAM2F code to generate an exit plane startline for both single and two phase 
motors which will interface with the JANNAF Standardized Plume Flowfleld 
(SPF) Code. 

The program is being used to generate nozzle and plume data for 
numerous applications such as: 

o Gas/ Gas-Particle Impingement (Heat-Transfer-Loads) 
o Rocket Nbzzle Performance (Thrust, I„_) 
o IR Signatures (Radiating Species) 
o RF Attenuation (Electron Densities) 

o Plume Radiation (Radiative Heat Transfer Gas/ Particles) 
o Vehicle Base Pressure 

o Base Heating (Convection-Recirculation), and 
o Flowfleld Properties for Contamination Predictions. 
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This three-volume report describes the program and Is Intended to 
provide the user the necessary information to use the code. It is 
anticipated that there may be some areas of the operation of the code that 
are not discussed In sufficient depth. Through feedback from the user the 
documentation will be periodicslly updated to reflect the experience of the 
various users. Additionally, coding errors which are uncovered will also be 
periodically provided to users. 

The code has been checked out and compared with a wide range of data. 
Validation of the program in the low density area of the plume has been 
verified with limited amounts of data. High altitude (vacuum) plume data 
are still required to further verify the applicability and range of use of 
the RAMP2F code. As these data become available the code will be further 
verified or modified to improve its applicability. 
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Appendix A 

SYMBOLS AND NOTATION 


Symbol 


Description 


f 

Kn 

M 

Nu 

Pr 

q 

r 

R 

Re 

T 

Greek 

r 

p 

p 

Subscripts 

C 

FM 

I 


drag coefficient, dimensionless 

incompressible drag coefficient defined by Eq. (A -21 

Knudsen number, dimensionless 

Mach number, dimensionless 

Nusselt number, dimensionless 

Prandtl number, dimensionless 

velocity, ft/sec 

particle radius, ft 

2 2 

gas constant, ft /sec /°R 
Reynolds number, dimensionless 
temperature, °R 

specific heat ratio, dimensionless 
viscosity, lbf-sec/ft^ 
density, slug/ft^ 

indicates the quantity pertains to t. = continuum flow regime 
indicates the quantity pertains to the free molecular flow regime 
indicates the quantity pertains to incompressible flow 


Superscript 

j indicates the quantity pertains to a particle species 


PRECEDING PAGE BLANK 


NOT FIIMED 
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APPENDIX A 

Particle drag and heat transfer quantities are calculated in terms 
of non-dimensional coefficients. These coefficients are functions of 
the local Reynolds number or the Knudsen number depending on the 
parameter being calculated. The following discussion summarizes the 
calculations employed in the gas-particle flow analyses. Details of 
the theory are given in the cited references. 

A.l PARTICLE DRAG COEFFICIENTS 

Particles in solid rocket motor nozzle and exhaust plume flow- 
fields encounter a wide range of flow regimes. Consequently any 
expression chosen for the drag coefficient must have the ability to 
take this into account. In a rocket nozzle and exhaust plume flow field 
the particles are continually being accelerated by the turbulent flow which 
has a temperature different than che particles. This situation is further 
complicated by rarefacJ.e; effects because the mean free path in the 
gas is comparable to dimensions of the particle boundary layer. 

Obviously the "standard drag curve" of a sphere as a function of 
Reynolds number (Ref. A-l) does not fit this description. This drag 
coefficient can only be applied to 

• Single solid sphere 

• Constant velocity 

• Still, isothermal, incompressible flow of infinite extent 

Li early two-phase studies the sphere drag coefficient was taken 
as the Stokes drag coefficient; 
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C 


D 


_ 24 
" Re 


(A>1) 


This expression, however, is valid only for Re L 1 if rarefaction 
effects are not taken into account. 

The Stokes drag coefficient is now used s a reference drag co- 
efficient. An incompressible drag parameter defined by the relation 


fJ * c D j '' c D J Stoke „: NP 


(A .2) 


was tabulated as a function of Reynolds number by Kliegel (Ref. A-2). 
This was corrected for rarefaction effects (Ref. A- 3) by the relation 


Cd _ c I (14 7.5 Kn) (1+2 Kn) + 1.91 K- 2 1 

Dj \ (1 + 7.5 Kn) (1+ 3 Kn) + (2.29+ 5.16 Kn) Kn 2 J 


(A. 3) 


The Knudsen number, Kn, (a measure of flow rarefaction) is related to 
the Mach and Reynolds numbers by 


KnJ = 1.26 \Jy~ mVR^, 


(A .4) 


M j = | Aq j |/V7RT • 


(A. 5) 


A more recent formulation for the drag coefficient is taken from the 
work of Crowe (Ref. A-4). The drag coefficient is defined as 
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C D = " C D I >^D + C D X 5 


where 


C D 


{ 


1.1 


1+1.1 (Kn) 




exp f 






1 - exp 


-Kn* 7 e- 


Kn(C 


D. 

me 

6 .6 



A.2 PARTICLE HEAT TRANSFER COEFFICIENTS 

At present there appear, to be no specific experimental heat trans- 
ier data for spheres of the particle sizes in the flow regimes encountered 
in a rocket nozzle- exhaust plume flow field. A review of semi-empirical 
expressions is given by Miller and Barrington (Ref. A-5). Until further 
work is done that is directly applicable to particle flows the continuum 
expression of Drake (Ref. A-6) 


Nujont = 2 + 0.459 (Re-*) 0,55 Pr 0,333 ; 


where 


Re^ = 2p | ATJ | rV p 

modified for rarefaction effects by Kawanau (Ref. A-6) 


Nu 


Nu = 


cont 


1.0 + 3.42 M/Re j Pr Nu 


cont 


appears to be the most applicable to two- phase flows. 
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Appendix 3 

SYMBOLS AND NOTATION 


Symbol 

h 

H 

P 

q 

8 

T 


Description 

local enthalpy, ft^/sec^ 
total enthalpy, ft^/sec^ 
pressure, lbf/ft^ 
velocity, ft/sec 
entropy, ft^/sec^/°R 
temperature, °R 


Greek 

q 

P 

A 

Subscripts 

o 

f 


Superscript 


weight flow ratio of oxidizer to fuel, dimensionless 
density, slug/ft 

denotes a change in a variable over a step length 

indicates the quantity pertains to the oxidizer 
indicates the quantity pertains to the fuel 


denotes a vector quantity 

denotes an average value over a step length 
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Appendix B 

The bulk of this document is concerned with the discussion of the 
analysis of gas-particle flows which have been shown to be non- isentropic 
as well as non-isoenergetic in the gas phase. However, one of the 
original objectives of the RAMP code development was to provide for a 
single phase (gas) calculation as an option. Computationally this 
entails simply considering the coefficients reflecting particle contri- 
butions to the compatibility equation as numerically equal to zero. As 
expected, there is no change in entropy and gas total enthalpy along 
gas streamlines, that is, the flow is seen to be isentropic and iso- 
energetic. However, for many liquid propellant motor applications a 
fuel rich mixture of the chamber combustion products is used to 
film cool the nozzle walls. This results in fuel striations across the 
nozzle where each streamline has a different energy level associated 
with it (i.e. non-isoenergetic flow). 

This case can be treated with a non -equilibrium chemical analysis 
by simply specifying the appropriate initial gas total conditions on each 
streamline. Recall that the pressure-density-velocity form of the 
compatibility equation is then used. Gas chemical concentrations are 
required for each streamline so that the gas thermodynamics are 
computed locally as the solution proceeds downstream in the nozzle- 
exhaust plume flowfield. However, if the flow is to be considered in 
chemical equilibrium or equilibrium/frozen, the applicable compatibility 
relations must reflect the non-isoenergetic nature of the striated flow 
condition. These relations are developed in the following paragraphs. 

To begin, the development of the species continuity equation 
(starting with equation 2*6) could be replaced by atomic conservation 
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equations. Moreover if the conservation of those atoms associated 
with the fuel and the oxidizer were considered then, for steady state. 


V * (p f q) = 0 (B*l) 

and 

V * (p 0 q) = 0 (B*2) 


would result. 

If the weight flow ratio of oxidizer to fuel (O/F ratio) is denoted by q 
then Eq. (B*l) and B»2) are satisfied if. 


q * Vq = 0 


(B-3) 


and 


V • (p$) = 0 (B.4) 

The assumptions inherent in arriving at Eq. (2*50) did not involve 
isoenergetic flow so that the momentum equation remains valid. The 
energy Eq. (2*105), however, remains valid only along each streamline 
and must be replaced by 

Jl 

h+f =Hfo) (B-5) 

The compatibility equations which apply along each Mach line are unaltered 
by the non-isoenergetic analysis since it is constructed based on the 
momentum equation and the global continuity equation (neither of which 
are altered by the non-isoenergetic flow phenomena). 
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The numerical solution to the governing equations is not greatly 
affected by the modification to the compatibility equation. In Eq. (3*146) 


A “l.2 


— (j^ + <* 1.2 **!.*) 
1.2 \ P 1.2 / 


while for a lower pressure boundary 


(s 3' s 1.2' *= 


1.2 \ P 1.2 / 


(B*6) 


(B*7) 


The finite difference analog to the streamline equation may be 
written: 


a. interior point 

n 3 = q 5 

b. an upper wall or upper free boundary point 

c. a lower wall or lower pressure boundary point 

q 3 =tl 2 

The oblique shock relations of Section 4 are unaffected since the 
discussion concerns a single streamline. All references to flow variables 
are understood to be a function of the O/F ratio. Since that ratio does 
not change across the shock the only alteration necessary is that of 
Eq. (4*11) which is altered to read 

P 2 ~ P ( s 2* q 2 ,T V’ P 2 = ^ s 2' q 2 ,1 'V t B * 8 ) 
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Since the Prandtl- Meyer calculation of Section 5 is also concerned 
with a single streamline, hence a constant value of i, it to is unaffected. 

It can readily be seen from the above discussion that the alterations 
necessary to incorporate the non-isoenergetic analysis are straight- 
forward from an analytic point of view. 
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Appendix C 

SYMBOLS AND NOTATION 


Symbols 


A 

B i 

c 

P 

H 

M 


N 

N 


P 

P 

q 

R 

S 

T 


u, v 
V 
w 
X 

Greek 


P 

0 


Description 

parameter defined by Eq. (2-46), 1/sec 
parameter defined by Eq. (2-103) 

2 2 o 

specific heat at constant pressure, ft /sec / R 
2 2 

total enthalpy, ft /sec 
molecular weight, lbm/lb mole 
number of moles 
number of discrete particles 
pressure, lbf/ft 2 
velocity, ft/ sec 

2 2 

gas constant, ft /sec / R 
entropy, ft 2 /sec 2 /°R 
temperature, °R 
velocity components, ft/sec 
velocity, ft/sec 
flow rate, lbm/sec 
position coordinate, ft 

density, slugs/ft 

inclination of the flow vector with respect to the x-axis, deg 
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SYMBOLS AND NOTATION (Cont'd) 


Subscripts 

c 

g.G 

m 


Superscript 


denotes chamber conditions 

indicates the quantity pertains to the gaseous species 

indicates the quantity pertains to the gas-particle 
mixture 

indicates the quantity pertains to the particle species 
indicates the quantity pertains to the particle species 
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Appendix C 

One oi the basic assumptions made in the development of equations 
describing the gas -particle flow is that there is no mass transfer between 
the phases. Consequently, calculations treating chemical reactions in gas- 
particle flows must reflect this assumption. Prozan (Ref. C-l) has shown 
for chemical equilibrium that the gas thermochemistry calculations can be 
uncoupled from gas-dynamic considerations. Gas thermodynamic properties 
(Ref. C-2) for chemical equilibrium or frozen flow are precalculated for nozzle 
chamber conditions. Tabulated thermodynamic data are then constructed in 
terms of the independent variables by expanding the flow from the chamber 
conditions. 

Thermochemistry calculations are, however, complicated by condensed 
species produced in the combustion process. Since the particles are inert, 
they perform no expansion work. Particle acceleration is accomplished by 
viscous drag forces exerted by the gas. This, coupled with the gas cooling 
resulting from expansion, produces dynamic and thermal lags. The work 
required to accelerate the particles combined with the heat exchanged between 
the particle and gas phases is a non-isentropic process as shown by the follow- 
ing relations! 


( V R) 

Tds ~ ■ * x 

Rpq cos 9 


NP 

dx = 0; 

j-1 


(C.l) 


dH - Tds + r 
P 





dx= 0. 


(C-2) 
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The term is a function of the velocity and thermal lags and in essence 
reflects the loss in work potential. Equation (C«l) indicates that if the 
algebraic sign on is positive then the change in the entropy level along 
a streamline will be positive. This is the case when the entropy in- 
crease resulting from the work required to accelerate the particle is 
greater than the heat transferred from the particles to the gas. Conse- 
quently there will be a corresponding decrease in gas total enthalpy. 

It should be noted that just the reverse is true if the increase in work 
potential due to heat transfer is greater than the work required to 
accelerate the particles. 

The original chemical equilibrium formulation (Ref. C-2) assumed 
that the gas and particles comprising the flow mixture were in dynamic 
and thermal equilibrium. Condensed species resisting from the com- 
bustion process were treated in the same fashion as the gas chemical 
species in that chemical reactions were permitted. The expansion pro- 
cess was then isentropic since gas-particle equilibrium was assumed. 
However, the correct thermochemical properties can be obtained if 
certain rr ations to the analysis are made. A step-by-step descrip- 

tion of th® modifications is contained in the following paragraphs. 

First the adiabatic flame calculation is made for the gas-part?., e 
mixture (Step 1, Fig. C-l). The mixture concentrations comprise the 
gas chemical species as well as the particulate species. Examination 
of the concentrations permits identification of the particulates. Since 
the assumption has been made that the gas and particulate phases 
do not exchange mass, the chemical reactions between the phases are 
suppressed by removing these from the list of possible reactions. 

Particle energy and mass are removed from the mixture (Step 2). The 
particle loading is then calculated from the relation 



M N 

W P-,,P 

M N 


wg 


( 03 ) 
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Input P. and 

Propellant 

Composition 

Step 1 




Fig.C-1 - Schematic of Two -Phase Thennochemical Calculations 
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Equation (C«3) gives the total particle mass loading relative to the gas. 
Individual particle loadings are obtained from empirical data for a given 
motor class and size. These data along with the particle enthalpy - 
temperature data (Ref. C-3) which defines the energy removed from the 
mixture become input data to the RAMP code. Next, the remaining gas 
energy and molar concentrations are then readjusted to obtain a one gram 
system (Step 3). The adiabatic flame calculation is made for the gas 
phase at the chamber pressure and temperature computed for the gas- 
particle mixture (Step 1). 

Equations (C« 1 ) and (C«2) indicated the gas phase expansion process 
to be non-isentropic . Hence the thermochemistry calculation must re- 
flect this condition. This is accomplished in the following manner. 

Equation (C-2) defines the change in gas total enthalpy as a function 
of entropy change and particle velocity lags. A schedule of gas total 
enthalpy changes can now be defined by the relation 

H G. = H c + AH i ; i= l ’ a (C*4) 

The values of AH^ are defined to cover the expected change in gas total 
enthalpy for a giver analysis. 

Figure C-2 shows a typical schematic of an H-S diagram for gas- 

particle flows. Each point Hq, Pc defines a particular point on the 

diagram from which the flow can be isen: pically expanding for the 

corresponding value of Sj. However, the flow is non-isentropic so that 

for each H/- , Pc an expansion from one or more states defined by a 
i 

change in entropy must be made. This is depicted schematically in 
Fig. C-2. Computationally this is accomplished by inputing to the 
code a schedule of AH^ and a corresponding change in entropy for each 
AH^. The chamber calculation is made for the gas phase (Step 3, 
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I 1 



= h c ±ah. 


| CEC | 



S 


Fig. C-2 - Schematic of Thermochemistry Table Construction 
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Fig. C-l) and an isentropic expansion from the reference conditions 
calculated for the number of entropy levels desired. These data are gene* 
rated in tabular form as a function of AH., ASL and Vfe respectively 
as indicated by the following schematic. This process is subsequently 


m i 



i 

i 

l 


AS. 

J 



Fig. C-3 - Schematic of Tabulated Thermochemistry Data 


repeated for each AHj. Since for a chemical equilibrium analysis the 
independent variables chosen are AH, AS and V the local thermodynamic 
and transport properties at a given point in the flow are obtained cy 
interpolating within the tables. 
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Appendix D 

SYMBOLS AND NOTATION 


Symbol 

Kn 

M 

N 

R 

Re 

Ro 

S 

T 

V 

V 

Greek 

y 

P 

P 

X 


Description 

Knudsen number, dimensionless 
Mach number, dimensionless 

number of collisions per molecule required for the gas 
to achieve relaxation 

gas constant, ft^/sec^/°R 

Reynolds number, dimensionless 

characteristic dimension of the flow field, ft 

distance along a streamline, ft 

temperature, °R 

gas velocity, ft/ sec 

mean molecular velocity, ft/sec 

specific heat ratio, dimensionless 

2 

viscosity, lbf-sec/ft 
density, slug/ft 
mean free path, ft 


PRECEDING PAGE BLANK NOT FILMED' 


D-iii 

LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-RREC TR 0867400-1 


Appendix D 

D.l DETERMINATION OF THE LOCAL KNUDSEN NUMBER, Kn 

The expansion of exhaust gases to very low back pressures results 
in a decrease in the gas pressure and density. This inturn produces a 
decrease in the number of intermolecular collisions. At some point, 
the effect of intermolecular collisions will cease to be a governing factor 
in determining the flowfield characteristics. The local Knudsen number 
(Kn) is used to indicate the local flow regime, i.e., continuum, vibrationally 
frozen, rotationally frozen or translationally frozen. The local Knudsen 
number is defined as the ratio of the mean free path (which is a measure 
of the average distance between intermolecular collisions) to a characte- 
ristic dimension of the flowfield. Kn is given by this relation (Ref. D-l): 

v 

In the relation for Kn, the mean free path is represented by "t?" A and 

the characteristic dimension is represented by T/-rjg-. The method 

employed here is to utilize the "sudden freeze" approximation (Ref. D-l) 

to determine when an energy mode freezes such that the gas temperature 

cannot relax fast enough to maintain the equilibrium rate of change. The 

RAMP program computes a xocal ".n at each mesh point and compares it 

til 

to l/Nj (the inverse of the collision number) to determine if the i 
energy mode has frozen. 

Assuming that the mean molecular speed V is determined by the 
Maxwell distribution function 
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V= (-|rt) 1/2 (Ref.D-2) 


were R is the particular gas constant. Introducing the local Mach 
number M, Kn can be written: 

Kn =(-|-r) / . dUaT)=^r • 

For rarefied gases the mean free path can be related to the viscosity \L 

X = (Ref.D-3) 

0.499pV 

— C 

where p is the local gas density. Let S = where R q is the character- 
istic dimension of the flowfield. Then dS = R q dS. With the above relation 
X = f(p), and introducing the characteristic dimension R^, Kn can be written: 


t m 2 / a \ 

Kn " 8 (.4995 \P VR 0 / 


But the local Reynolds number Re is: 


d(fn T) 

d3 


Therefore, 


Kn = 0.25050 ir 


d(in T) 

d3 


As a first approximation to the derivative » let 

d5 
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d(ln T) 
dS 


to (Ti/ 


(Tl/T 2 ) 

o 


so that the local Knudsen number, Kn as coded is given by: 


Kn = 0.25050 n 

The characteristic mesh is shown below to illustrate the flowfield 
variables used to determine the local Knudsen number 



D.2 NON-CONTINUUM FLOW REGIME CRITERIA 

The criteria used by the RAMP program to determine the "freezing" 
of particular energy modes is to compare the local Knudsen number (Kn) 
with the reciprocal of the collision number. When Kn is equal to or 
greater than the reciprocal of the collision number for the vibrational, 
rotational, or translational modes that particular mode is considered to 
be frozen. The freezing of the first two energy modes (vibrational and 
rotational) represents a transitional computational region from the con- 
tinuum to free molecular regimes. During this transitional period the 
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gaseous flowfield calculation is performed using the characteristic 
equations with the thermodynamic properties modified to reflect the 
respective flow regime. That is, the gaseous ratio of specific heats 
is altered to reflect the transition to a non-continuum regime. Once 
the translational energy mode has frozen the flowfield becomes "free 
molecular." In this case the solution becomes a source flow solution 
with a constant gaseous ratio of specific heats. 

D.2.1 Non-Continuum Flowfield Gas Properties 

The energy modes of the gaseous flowfield freeze one at a time 
in order of increasing energy requirements. The vibrational mode 
freezes first, the ratio of specific heats (gamma) is set to 7/ 5 and the 
flowfield solution proceeds in the conventional streamline normal fashion. 
See Ref. D-4 for the characteristics equations. The rotational mode 
freezes next in which case gamma is computed as a function of 
Kn (third order polynomial) and varies from 7/5 to 5/3. This flowfield 
solution is not altered except for the determination of gamma. The 
last energy mode to freeze is the translational mode. Once the trans- 
lational mode has frozen, the flowfield is considered to be in the free 
molecular regime. Gamma is set to 5/3 and is held constant. The 
effect of the freezing of the three gaseous internal energy modes on the 
gas ratio of specific heats is depicted pictorially in Fig. D-l. 

Once the flowfield calculation has switched to the free molecular 
regime the gas streamlines remain straight and the gas velocity which 
is directed along streamlines remains constant. The gas constant, 
temperature, and flow angle also remain constant. The gas density 
varies inversely as the streamtube area and the gas pressure is 
obtained from the equation of state. 
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Local Knudsen Number, Kn 

Fig.D-1 - Variation of Gamma with the Local Knudsen Number 


D.2.2 Non-Continuum Flowfield Particle Properties 

If particles are present in the free molecular flowfield their status 
is computationally very similar to the gas. Particle velocity, tempera- 
ture, and enthalpy remain constant. The particle flow angle remains 
constant, and the particles are directed along straight streamlines. The 
particle and gas streamlines are not necessarily parallel. The particle 
density varies inversely as particle stream tube cross-sectional area. 
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Appendix E 

In flow problems where the gas may be considered in equilibrium 
(chemical and thermodynamic) at every point, two parameters are sufficient 
to define any of the other thermodynamic variables, either by assuming a per- 
fect gas or by utilizing the results of chemical equilibrium calculations of the 
gases involv *d. If the gases cannot be considered to be in equilibrium, a 
r mt-by-poi it evaluation of its composition via int< ^rating the individual 
species continuity equations is required. This appendix addresses the 
procedure utilized in the RAMP program to perform such calculations. 

A detailed description of the rate processes that occur in rocket ex- 
haust P-.V.-5 requires that a myriad of mechanisms he considered to include 
all the possible chemical reactions of dissociation, formation, recombination, 
etc. 


All of these, however, can be treated with a very general formalism. 
In the form usually quoted in chemical kincii.cs (Ref. E-l) the phenomeno- 
logical law of mass action states that the rate of a reaction is proportional 
to the product of the concentrations of the reactants. Thus, for a general 
reaction of the form 



(E.ll 


the net rate of production for any participating species for which the 
stoichiometric coefficients f/! and |/|' ar*- not equal cap then be written as 


w = k 


N V\ N „ 

f /7(C ) l -k, |7 (C.^i 
i=l D 


(E.2) 
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where C. is the species concentration defined as = p F. with p being the 
density and being the species mol/mass ratio, respectively. Assuming 
small deviations from equilibrium, the forward and backward reaction rate 
constants, kj and k^, respectively, can be related to the concentration equi- 
librium constant and to the pressure equilibrium constant as follows: 


N 




K ( 3* T) 

Jr 


i=l 


ivl 


V'P 


(E.3) 


The significance of the pressure equilibrium constant K is that it can be 

F 

easily evaluated for any reaction using tabulated values of Kj the equilibrium 
constant for formation from the elements. Values of Kj are commonly tabu- 
lated in conjunction with specific heats, entropies and enthalpies as a function 
of temperature, and are available in general for most species. An equally 

convenient method exists for determining K from the change of free energy 

P 

accompanying the reaction, i.e.. 


K p = expf-AG/KT) 


(E.4) 


where AG is the change in free energy during the reaction process. Free 
energy values are also available for most species in tabular form. This 
method is used to compute K p in the RAMP program. 


Using Eq. (E.3) to eliminate the backward rate constant from Eq. (E.2) 
the general production rate equation can be written as 

N 

W - vi) 


w = k. 


N 

n (c { ) 


n 


LIT) 


i = l 


K 


N 

77 (C.) 

» « ^ 

1 = 1 


vi' 


(E.5) 
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Finaily, the net rate of production for any species participating in the 
reaction, either as a reactant or as a reaction product, is then given by 

Wj = (V[' - u{) w (E.6) 

Since most reaction systems involve numerous simultaneous reactions, 
the net rate of production of species i usually equals a sum of terms. Thus, 
for an aribitrary number of M simultaneous reactions, all involving species i, 

M 

w. = ^ k k = l....,M (E.7) 

k=l 

where w. . denotes the net rate of production of species i due to reaction k 

1, K. 

only. 


For reasons of computational speed and efficiency, the program con- 
tains explicit expressions, as obtained from Eq. (E.5), for the most commonly 
encountered reaction mechanisms. 1 /elve types of reaction mechanisms are 
considered as possible contributors to the calculation of the net rate of pro- 
duction, w.: 

Reaction type 


(1.7) 

A + B 

— * 

C +D 

(2, 8) 

A + B + M 


C + M 

(3,9) 

A+B 


C + D + E 

(4, 10) 

A+B 


C 

(5.11) 

A + M 


C + D + M 

(6. 12) 

A + M 


C + M 


Reaction types (7) through {12} Correspond to reaction types (i) through 
(6), but proceed in the forward direction only. 

To reduce roundoff and truncation errors, the forward and backward 
rates for each reaction are computed separately. All contributions to the 
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molar rate of production of a given species are then computed and added 
algebraically to form matrix coefficients (discussed later). Since reaction 
types (7) through (14) proceed in the forward direction only, the second 
term on the right-hand side of Eq. (E.5) is disregarded in calculating the 
contributions to the coefficient matrix. 


In reactions (2), (5) and (6) as well as in (8), (11) and (12), M denotes 
a third body species which can be specified. For these reactions the situation 
often occurs where for various third bodies the respective rate constants 
differ only by a constant multiplier. These multipliers can be considered as 
third body efficiencies or weighting factors. If such a case is encountered, 
the third body species mole mass ratio becomes effectively a fictitious 
mole mass ratio, consisting of the weighted sum over all those species having 
a nonzero weighting factor, i.e., 

F M ■ E‘i F M. < E * 

l l 


where f. are the weighting factors. 

The forward rate constant, k^, is generally a function of temperature. 

It is most often expressed in Arrhenius form. Again, for speed and efficiency 
in computation, the rate constants are divided into five types: 

Rate Constant Type 


= A exp(B/ SI T) 

= AT~ N exp(B/3T) 

= AT' N exp(B/$T M ) 


The equations presented in this appendix provide a very general forma- 
lism for the evaluation of various rate processes. The specification of par- 
ticular systems and associated rate constants is up to the program user. 
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Cor sidering now the general species continuity equation 

pq. VC. = w (E.ll) 

and making use of the foregoing discussion of the rate process we now proceed 
to de;ciibe a calculational technique for determining the individual species 
composition m a point - by -point basis. The description of this process is 
substantially simplified if Eq. (E.5) is specialized to a particular reaction 
type, say lumber (7) from Eq. (E.8) which is a one-way, two-body reaction. 

A + B — ^ C •*- D (E.12) 

the net production rate for this process is 

w = k f p 2 F A F B (E.13) 

auu the species continuity equation for species B then becomes 

#.VF b = k f p 2 F A F B (E. 14) 

which along streamlines becomes 

0F B 2 

w-ar = f a f b (E15 > 

This equation can readily be solved using finite difference techniques 

employing explicit relationships, such as Euler or more sophisticated schemes, 

such as Runge-Kutta. The step size for integrating this equation however is 

severely limited by stability criteria. It can be seen from Eq. (E.15) that the 

rate of change of a species along the streamline becomes increasingly larger 

as the flow speed is slowed, too density increased, or for fast reaction rates. 

In rocket engine problems, combinations of slow speeds, high densities and 

fast reaction rales (i.e., quasi -equilibrium) are quite common and integration 

-8 

step sizes so small (i.e., < 10 meters) are encountered that the solution 
becomes impractical in terms of computation time. 
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For this reason, the technique described in Ref. E-2 based on a lineari- 
zation of the production rates was utilized. Writing Eq. (E.15) in finite differ- 
ence form over a streamline step from station n to n+1. 


F 


n+1 



(E.16) 


And evait. all the species concentrations at the downstream point results 

in a set of simultaneous nonlinear algebraic equations. In order to solve 

these equations we must then linearize the term F . F_ which is accom- 

n+1 n+1 

plished following the lead of Ref. E-2. If this term is expanded in terms of 
its values at station n along with the increments over n to n-i-1 we can obtain 
the following expression. 





+ F B 



(E.17) 


neglecting products ot differentials which are assumed to be of second order 
importance. Equation (E.16) can now be written in its linearized form. 

Let C = AS p/q 



F b „ 

= f b 

- c 

f a f 


n+1 

n 


n 

and 

F a 

n+1 

= f a 

n 

- c 



+ F F 

B B A _ . 

n+1 n+1 


F B + F B . F A 
n+1 


(E. 18) 


Equation (E.18) can then be expressed in terms of a set of unknowns and 
calculable coefficients, C. Rewriting these we obtain 


n+1 


n+1 


f a - cf a < f b_J- cf b < F A > 


n n+1 


n “n+1 


(E.19) 


f b - C F B < F A > ‘ C F A ^ F B > 
n n n+1 n n+1 
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F a dtCF B ) t (CF A )F B ^ = F a 


n+1 


n n+1 


f A (CF b )M1-CF a) F b = F 
n+1 B n A n n+1 n 

A matrix can now be formed using totally known information. 


1 +CF b 

cf a 


" F A 


" f a" 

n 

n 


n+1 


n 

cf e 

n 

1 - cf a. 


f b 

_ n+l_ 


F B 

- 


(E.19) 

Cont'd 


(E.20) 


The matrix [A] [x] = [b 3 is then solved for the unknown compositions F a , 

n+1 

F_ via a triangulation technique. Although consuming more time per 
n+1 

integration step than an explicit formulation, the implicit technique employed 
here is unconditionally stable permitting much larger step sizes, thus allowing 
solutions to be obtained for problems where the small steps required by the 
explicit technique prevented even the consideration of the case. Finally it 
should be recalled that an extremely simple case was chosen only for purposes 
of illustration and the general technique coded in RAMP handles many species 
(see input guide) with multiple reactions. 
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